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■ Strongly correlated Fermi system plays a fundamental role in very different areas of physics, 

from neutron stars, quark-gluon plasmas, to high temperature superconductors. Despite the broad 
| applicability, it is notoriously difficult to be understood theoretically because of the absence of a 

£Sj . small interaction parameter. Recent achievements of ultracold trapped Fermi atoms near a Feshbach 

resonance have ushered in enormous changes. The unprecedented control of interaction, geometry 
O ' and purity in these novel systems has led to many exciting experimental results, which are to be 

urgently understood at both low and finite temperatures. Here we review the latest developments 
of virial expansion for a strongly correlated Fermi gas and their applications on ultracold trapped 
Fermi atoms. We show remarkable, quantitative agreements between virial predictions and various 
recent experimental measurements at about the Fermi degenerate temperature. For equation of 
'■ state, we discuss a practical way of determining high-order virial coefficients and use it to calculate 

accurately the long-sought third-order virial coefficient, which is now verified firmly in experiments 
at ENS and MIT. We discuss also virial expansion of a new many-body paramter - Tan's contact. 
We then turn to less widely discussed issues of dynamical properties. For dynamic structure factor, 
£^ ■ the virial prediction agrees well with the measurement at the Swinburne University of Technology. 

For single-particle spectral function, we show that the expansion up to the second order accounts 
for the main feature of momentum-resolved rf-spectroscopy for a resonantly interacting Fermi gas, 
as recently reported by JILA. In the near future, more practical applications with virial expansion 
1 are possible, owing to the ever-growing power in computation. 
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I. INTRODUCTION 



Universal strongly correlated Fermi systems: From dilute neutron matter to ultracold trapped Fermi 

atoms 



The strongly correlated Fermi gas is a ubiquitous system in nature It appears in the quark-gluon plasmas in the 
early Universe [2], neutron stars 0,01, high-temperature superconductors and most recently in ultracold atoms 
[6H8| (see Fig. QJ. The strong correlation is a result of a large separation of length scales and the Fermi system is close 
to an interesting universal limit with infinitely large scattering length and zero effective range of interaction 0, [To| . 
The absence of length scale implies that the type and detail of interactions are not important. It is anticipated that 
universal behaviors in both static and dynamic properties would emerge [loMT2j . 
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Figure 1: (color online) Ubiquitous strongly correlated Fermi gases in nature. 

Dilute neutron matter is a good example of strongly-correlated Fermi systems [H, L3| • The neutron scattering length 
is about a s ~ —18 fm and the effective range is r$ ~ 2.8 fm <C a s . For typical neutron densities 0.1/9/v > P > 10 _4 /3at, 
where ~ 0.16 fm~ 3 is the saturation density of nuclear matter, the dimensionless interaction parameter Uf \a s \ 3> 1 

1/3 

while kp \ro\ is small. Here, kp — (3ir 2 p) is the Fermi wave-vector. Therefore, the neutron matter is close 
to the unitary limit, with which the s-wave scattering amplitude becomes saturated at a zero-energy resonance. 
Understanding the nuclear matter appears to be a challenging many-body theoretical problem fl3l [l4| . 

In this context, a unitary atomic Fermi gas realized recently in ultracold atom laboratory attracts particular atten- 
tion It serves as a new paradigm for studying strong-correlations because of its unprecedented controllability 
and purity. By tuning an external magnetic field across a collisional Feshbach resonance [l5(, the interatomic attrac- 
tions in a two-component Fermi gas can be changed precisely from weak to infinitely strong, leading to the observation 
of crossover from Bardeen-Cooper-Schrieffer (BCS) superfluids to Bose- Einstein condensates (BEC) [T|| [Tt]], which 
was anticipated long time ago |l8M20| . At the resonance, the s-wave scattering length a s is exactly infinity and the 
effective range of interaction is negligible. This unitary limit can now be routinely achieved in laboratories with 
fermionic potassium-40 ( 40 K) [lT| and lithium-6 ( 6 Li) atoms [l7| . 

Experimentally, the near resonance regime was first approached by O'Hara et al. in 2002 with 6 Li atoms [2l| . 
The stability of atomic Fermi gases under strong attractions was observed and the ground state energy was found 
to reduce significantly with respect to that of an ideal, non-interacting Fermi gas. Since then, a number of different 
aspects of a unitary Fermi gas have been characterized after substantial experimental efforts. Hydrodynamic ballistic 
expansion and collective excitations due to strong correlations were confirmed (22Tt2^ |: superfluidity at BEC-BCS 
crossover was unambiguously verifie d by ge nerating quantized vortices (2oT | ; universal thermodynamics was evidenced 
by heat capacity measurement (Til |26h321| ; nearly ideal hydrodynamic flow and universal viscosity was observed 
[33|. Recently, these studies have been extended to Fermi gases with unequal densities for the spin- up and spin- 
down components [U, [35| and to Fermi gases in low-dimensions [36T- |38| | . giving the prospects of realizing exotic 
inhomogeneous Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superfluidity [39r - l42| and Berezinskii-Kosterlitz-Thouless 
(BKT) transition [43T - [4o1 |. Along with rapid experimental progress, new measurement techniques have been developed 
to characterize strongly correlated Fermi gases. These include the momentum-resolved rf-spectroscopy for measuring 
single-particle spectral function (Trl l47j and the two-photon Bragg spectroscopy for dynamic and static structure 
factors [H|. 

In contrast, the parallel theoretical development is much slower. There are numerous activities on developing better 
strong-coupling theories (2(1 149| - |58| or ab-initio quantum Monte Carlo (QMC) methods (oTSWriil ]. However, a deep 
understanding of strongly correlated Fermi gases is prohibited because of the absence of a controllable small interaction 
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parameter. The use of standard strong-coupling theories requires infinite order expansions and the truncation to a 
particular order can not be fully justified a priori [6^,[6^|. At this stage, numerically exact QMC simulations are less 
accurate than one might expect, suffering from either the notorious sign problem for fermions [6l| or the finite size 
effects in small samples used in the simulation (60l. |62|. 

In this respect, exact results of strongly correlated Fermi gases in some non-trivial limits are very valuable. Two 
recent efforts are notable. In the limit of short-distance, large- momentum, and/or large- frequency, Tan derived a 
set of exact universal relations (67l [68[ . It was shown that all the limiting behaviors are governed by a many-body 
parameter called the contact, which measures the density of pairs within short distances. Tan's relations can be 
conveniently understood using the short-distance and/or short-time operator product expansion (OPE) method (69j . 
which separates in a natural way the few-body physics from many-body physics. In another limit of high temperature, 
quantum viriai expansion [70l - l79| provides another rigorous means to bridge few-body and many-body physics. The 
properties of a strongly correlated Fermi gas, either static (70l. l72l l73l 1771 r78| or dynamic [75l. [7a|. can be expanded 
non-perturbativel y using some exact expansion coefficients or expansion functions, which are calculable from few- 
fermion solutions [80IH86| . Both Tan relations and viriai expansion give useful insights into the challenging many-body 
problem, though in the different perspective. 

In this paper, we review the recent theoretical development on quantum viriai expansion, and show that viriai 
expansion gives a complete solution of strongly-correlated Fermi gas above the Fermi degenerate temperature. We 
focus our attention on ultracold atomic Fermi gases, and compare in a quantitative way the viriai expansion predictions 
with available experimental measurements for various fundamental properties. We note that the viriai expansion has 
also been used frequently to study the equation of state of neutron matter [87T - |9^ |. 



B. Overview of viriai expansion picture 



Quantum viriai expansion, alternatively referred to as quantum cluster expansion, is a standard method in quantum 
statistical mechanics j93l. |94(|. It is practically useful for a dilute quantum gas. The basic idea of viriai expansion is 
simple. Though we have a strongly correlated system at low temperatures, with increasing temperature the correlation 
between particles would become increasingly weak. At sufficiently high temperatures, the scattering cross section is of 
the order the square of the thermal de Broglie wavelength, which becomes much smaller than the average inter-atomic 
distance. As a result, the inclusion of few-body correlations is already sufficient to describe the underlying properties 
of the system. These few-body correlations can be exactly taken into account using few-particle solutions and viriai 
expansion. 

As a concrete example, let us consider the thermodynamic potential Q for a given Hainiltonian "H, which in the 
grand canonical ensemble is given by (95| . 

fl = -k B TlnZ, (1) 

where k B is the Boltzmann constant, 

Z = Tr cxp [- (H - MO /k B T] (2) 

is the grand partition function, and M is the field operator of total number of particles. The thermodynamic potential 
can be written in terms of the partition function of clusters, 

Q n = Tr„ [cxp (-UJk B T)] , (3) 

where the integer n denotes the number of particles in the cluster and the trace Tr„ is taken over n-particle states 
with a proper symmetry. Q n is calculable using the complete solutions of a n-particle system. The grand partition 
function then takes the form 

Z = I + zQi + z 2 Q 2 + z 3 Q 3 • • • , (4) 

where z = cxp(/i/k B T) is the fugacity (95| . At large temperatures, it is well-known that the chemical potential \x 
diverges to — oo, so the fugacity would be very small, zCl. By Taylor-expanding In Z in powers of the small fugacity, 
it is obvious that 

n = zfi« + z 2 n< 2 > + + • • • , (5) 

where can be expressed in terms of Qi (i < n) and therefore contains the contribution from few-body physics up 
to n-particles. 
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In principle, all the properties of a quantum gas could be cluster expanded in powers of fugacity, no matter how 
strong the interactions. The fugacity is a natural small parameter at large temperatures. Naively, virial expansion 
is applicable when z < 1. For a two-component spin- 1/2 Fermi gas, using the fact that the fugacity is roughly 
equal to the phase-space density pA^ B /2, where p is the density, \ c ib = [2,7rh 2 / (m^T)] 1 / 2 is the thermal de Broglie 
wavelength and m the mass of atoms, one can estimate that virial expansion is useful at temperature T > Tp. Here 
Tf = H 2 k F / (2mkB) is the Fermi degenerate temperature. 



Q( 2 ) = 





u + 




Figure 2: (color online) Diagrammatic representation of the contribution of two- particle scattering process to the thermodynamic 
potential. Here Ti is the two-particle vertex function, obtained by summing all the ladder-type diagrams. The dashed line and 
solid line represent the bare contact interaction Uo and the single-particle Green function, respectively. For details, see refs. 
Ei and [H. 



Though virial expansion is a large-temperature expansion, it has an intrinsic relation with low-temperature strong- 
coupling diagrammatic theory. In the absence of a small interaction parameter, we may conjecture that a reliable 
strong-coupling theory of strongly correlated systems should be developed by successively including few-particle 
scattering process. Thus, we may write 



Q = Q« + + n (3) + 



(6) 



where is the thermodynamic potential of a non-interacting system, and fiW (n > 2) is the contribution from the 
n-particle scattering process, which is to be calculated at all temperatures by summing a series of diagrams to infinite 
order (i.e., the n-particle vertex function T n ). We shall refer to such an expansion as the diagrammatic expansion. As 
an example, in Fig. [5] we show the dia gram matic representation of Q( 2 \ It sums up all the two-particle scatterings 
via the two-particle vertex function [18j]. In the language of functional path- integral method, f^ 2 ' corresponds to 
the gaussian fluctuations around the mean- field saddle point [I8I - F20L |53| . It is obvious that at large temperatures, by 
expanding fiw (i < n) in powers of z, we can calculate directly fl^ appeared in the virial expansion. In this respect, 
virial expansion and diagrammatic expansion are closely related. Both of them are the expansion in few-particle 
correlations. The advantage of the diagrammatic expansion is that it takes into account the few-particle scatterings 
in the medium, and therefore is applicable at all temperatures. It is a natural generalization of virial expansion to 
the low-temperature regime. These two expansion theories are sketched in Fig. [3J 

Ideally, for a strongly correlated system, we anticipate that f^ n ) becomes less important with increasing n and the 
diagrammatic expansion thus converges. Indeed, for a unitary Fermi gas at the BEC-BCS crossover, the theoretical 
calculation of Sl^ at zero temperature gives fairly accurate equation of state [53[, as confirmed by the latest experi- 
mental measurement [§6j|. Others terms of fi( n ) with n > 3 are notoriously difficult to obtain, but are conjectured to 
be small at low temperatures. Virial expansion provides systematic determinations of fl( n > at high temperatures and 
may shed light on their low temperature behavior. 

At this point, we may anticipate that the applicability of virial expansion is not limited to small fugacity z < 1. 
The expansion could be meaningful in the deep quantum degenerate regime through an analytic continuation across 
the point z = 1, and therefore is applicable down to the superfluid phase transition temperature T c . The pursuit of 
such an improved virial expansion theory is a theoretical challenge. 
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Figure 3: (color online) Schematic illustration of virial expansion and diagrammatic expansion, both of which are expansion in 
few-particle correlations. It is desirable to find an improved theory, which connects smoothly these two expansion theories. 



C. Key technical issues in the latest development of virial expansion 

Despite the usefulness of virial expansion, its application to strongly correlated quantum gases is less documented 
in the literature. There are two severe technical difficulties in applying virial expansion to a homogeneous system: (i) 
insufficient knowledge on the exact few-particle solutions and (ii) continuous energy spectrum. As a result, it seems 
impossible to calculate the essential few-particle cluster partition function Q n when n > 3, as the calculation requires 
infinitely large number of energy levels. Therefore, previous applications of virial expansion have been restricted to 
the second order, where Q2 can be calculated using an elegant phase-shift formalism derived by Beth and Uhlenbeck 
in 1937 [93, H. 

The latest development of virial expansion, to be reviewed in this paper, relies on the recent theoretical progress 
on the exact few-particle solutions of trapped strongly interacting fermions [8cM86j . Due to the trapping potential, 
the energy spectrum becomes discrete. As the thermal energy fc^T provides a natural high-energy scale in the cluster 
partition function, the number of energy levels required by the calculation is finite. In principle, we can always 
calculate numerically these energy levels using the ever-growing computation power, if few-particle solutions are not 
known analytically. In addition, virial expansion of dynamic properties becomes possible, based on the calculated 
few-particle wave- functions (75l . [7o| . 

To get back to the homogeneous system, one can utilize the so-called local density approximation, which treats 
the trapped Fermi gas as a collection of many locally uniform blocks. In the unitary limit, it is found that the virial 
expansion results for trapped and homogeneous systems are convertible by some universal relations [13 , [lH . The 
details will be discussed later. 



D. Model Hamiltonian 



Throughout this Review, we focus on the strongly-interacting Fermi gases with zero-range interactions in three 
dimensions, which have emerged as the simplest strongly-correlated model system that has been accessed experimen- 
tally with ultracold atoms of 6 Li and 40 K [g-Q- They also serve as a "bare bone" description of nuclear and neutron 
matter. The generalization of virial expansion in low dimensional systems is straightforward [74j . 

The use of zero-range interactions can be understood from Fig. 01 which plots the short-range behavior of the 
zero-energy scattering wave- function i/i re i(r) for real interatomic interactions with a range of interactions ro- For an 
ultracold dilute Fermi gas, ro is typically at the order of 10 -9 m, much smaller than the mean inter-particle distance 
p -1 / 3 ~ 10~ 6 m. As a result, the complicated short-range behavior of the wave-function, which corresponds to 
high-energy physics, becomes irrelevant for typical atomic collisions. Therefore, we may set effectively tq = and 
approximate ip re i (r) oc l/i — 1 /a s at r ~ ro = [§4j . Here a s is the s-wave scattering length [94| . This is the so-called 
Bethe-Peierls (BP) boundary condition, which is equivalent to the s-wave zero-range interactions (or the so-called 
pseudopotential j94(). We note that, when the s-wave scattering length becomes positive, the interaction potential 
will support a two-body bound state with energy Eb = —h 2 /(ma 2 ). Therefore, in the unitary limit, where a s — > ±00, 
a shallow two-body bound state with infinitely small energy emerges. For more details, see ref. (la |. 

In the ultracold atom experiments, a harmonic trap is necessary to prevent atoms from escaping. We thus consider 
N fermions in a three-dimensional isotropic harmonic trap Vt(x) = mu:\{x 2 + y 2 + z 2 )/2 with the same mass m and 
trapping frequency cjy, occupying two different hyperfine states or two spin states. The zero-range s-wave interaction 
between fermions with unlike spins is replaced by the BP boundary condition. That is, when any particles i and j 
with unlike spins are close to each other, = |x, — Xj| — > 0, the many-body wave function ip (xi,X2, ...,xjv) with 
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Figure 4: (color online) Explanation on the use of zero-range interactions. For a dilute quantum gas at low temperatures, the 
short-range behavior of zero-energy scattering wave-function between spin-up and spin-down fermions (represented by different 
colors) can not be accessed. It can therefore be safely replaced by a simple asymptotic form, ri[> re i{r) oc l — r/a s , which defines 
an s-wave scattering length, a 3 . The interatomic interactions are characterized by the single parameter a a only. 



proper symmetry should satisfy |99j, Il0fj | , 

ip = Aij(X-ij 



i{Xfe^i,j}) 



1 1 



(7) 



where _/4.y(Xy, {x^j}) is a function independent of r^. This BP boundary condition can be equivalently written 



lim 



r M - >-o dr i:j a s 
Otherwise, the wave function ip obeys a non-interacting Schrodinger equation, 

* r *a 



E 



K 1 1 

w2 i 2 i 2 i 2 i 2\ 

2m 2 



ip = E-iJj. 



(8) 



(9) 



In the second quantization, the system can be instead described by the model Hamiltonian, 



u = E 

<T=t4 



dx^(x) 



2m 



+ y T (x) - Ma 



'0a(x) + C/o / dxipl(x)xjjUx)ijji(x)tp^(x), 



(10) 



where the chemical potentials ^ and /ij, could be different due to unequal spin-populations. The zero-range interaction 
is given by a contact potential Uq5(x — x'). The bare interaction strength Uo is to be renormalized by two-particle 
vertex function T 2 in the vacuum [95|], using 



1 



1 E 



(11) 



where the momentum k has a high-energy cut-off, k < A = r$ , in accord with the use of zero-range interactions, fTo 
scales to zero as the cut-off momentum A — > oo. 

The (single-channel) model Hamiltonian shown above provide the simplest description of ultracold 6 Li and 40 K 
atoms near broad Feshbach resonances [6-8]. In case of narr ow Feshb ach resonances, it is necessary to use a two- 
channel model and to include molecules in the closed channel [lOll . Il02| . 
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E. Brief introduction to Tan relations 



Here we introduce briefly the exact Tan relations, which were derived by Shina Tan in 2005 [63, [fig. These Tan 
relations link the asymptotic behavior of many-body systems at short-range, large-momentum, and high-frequency to 
their thermodynamic properties. For instance, the momentum distribution p rT (q) falls off as I/q 4 at large momentum, 
the pair correlation function fffj^Xj — Xj) = J dX-ij (pf(xi)p±(xj)) diverges like 



9n( r 



0) 



16tt 2 



a,r. 



(12) 



S' IJ 



and the rf-spectroscopy has the tail of I/u> 5 ^ 2 at large frequency |l03l . Il04l |. All the Tan relations are related to 
each other by a single coefficient X, referred to as the integrated contact density or contact. The contact measures 
the probability of two fermions with unlike spins being close together [6j|. It also links the short-range behavior to 
thermodynamics via the adiabatic relation, 



dE 



|0(-l/a 8 ) 



S,N 



Airm 



(13) 



which gives the change in the total energy E due to adiabatic changes in the scattering length. The fundamental 
importance of the Tan relations arises from their wide applicability. They are useful at both zero or finite temperature, 
superfluid or normal phase, homogeneous or trapped, few-body or many-body systems. 

While the original rigorous derivation by Shina Tan is difficult to follow, the underlying physics of Tan relations 
can be easily understood from several points of view (68| . The simplest way is from the two-body wave function under 
the BP boundary condition, tp re i(r) oc 1/r - l/a s . Naively, the momentum distribution p a {q) is simply the square 
of the Fourier transform of ^> r ez( r ) an d the pair correlation function g^\(r) oc |^/> re ;(r)| . The asymptotic behavior of 
p a (q) oc q and g^{r) oc r~ 2 is then straightforward to check. 

At the many-body level, Tan's relations can be elegantly proved by using the short-distance and/or short-time 
operator product expansion (OPE) method (d9| . in which the contact I is identified as 



1 = XJl / dxV4(x)V>j(x)^(x)</> t (x). 



For example, the adiabatic relation Eq. 
model Hamiltonian 1691, 



(14) 



l3| can be obtained directly by applying Hellmann-Feynman theorem to the 



dE 



d(-l/a s ) 



S,N 



on 



d(-l/a s )/ d{-l/a s ) 



Atttyi 



(15) 



The last step follows the renormalization for the bare interaction strength Eq. (fTTj) . 

The contact is a fundamental parameter that characterizes the many-bo dy properties of strongly correlated Fermi 
gases. Recently, its measurement receives considerable attentions |l05l4l08t. It turns out that the most accurate way 
is through the Tan relation for spin-antiparallel static structure factor |109| . which is obtained by a direct Fourier 
transform of pair correlation function, 



Sn (l > k F) 



1 
Wq~ 



1 



na s q 



(16) 



The simple power- law tail of 1/q in the structure factor relation is more amenable for experimental measurement than 
the q~ 4 or w~ 5 / 2 tail in the momentum distribution or in rf-spectroscopy. In the latter two cases, the fast decay due 
to the higher-order power law imposes more stringent signal-to-noise requirements at a given momen tum or frequency. 
Experimentally the static structure factor can be measured by two-photon Bragg spectroscopy [H, Il06| . 



F. Brief summary of virial expansion results 



We now summarize briefly the main results of the latest development in virial expansion. In general, thermodynamic 
properties such as the thermodynamic potential can be expanded in terms of some virial coefficients, while dynamic 
properties, i.e., the single-particle spectral function and dynamic structure factor, can be expanded in terms of some 
virial expansion functions. The latest developments of virial expansion include: 
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(i) The third virial coefficient 63 for thermodynamic potential has been precisely determined [72l . l78| . For a homo- 
geneous Fermi gas in the unitary limit, it is found that A63 = 63 — = —0.3551030264897 [78[. Here b^p is the third 
virial coefficient of an ideal, non-interacting Fermi gas. This theoretical prediction has been confirmed experimentally 
with the expe rimental value Ab3, e xtp = —0.35 ± 0.02 [3Cj and has been independently checked by a field theoretic 
approach IllOl . which gives A&3 = —0.3551 ± 0.0001. The fourth virial coefficient in the unitary limit has also been 
calculated [78] . but with much less accuracy, A64 = — 0.016 ± 0.004. The virial coefficients of a trapped system b n .T 
and a homogeneous system b n are related by [z2, & rl ,T = n~ 3 / 2 b„ (n = 1, 2,3, ■ • ■ ). These virial coefficients predict 
an accurate equation of state for a trapped Fermi gas in the unitary limit, for temperature down to 0.5TV [66| . We 
review the virial expansion of thermodynamics in Sec. II. 

(ii) The contact parameter I can be virial expanded in terms of contact coefficients c n 1771 . F or a homogeneous 
Fermi gas in the unitary limit, it is predicted that c 2 = 1/tt and c 3 = —0.1399 ± 0.0001 [11C( . In analogy with 
the virial coefficient, the contact coefficients of a trapped system c n ,T and a homogeneous system c„ are related by, 
c n ,T = n~ 3 / 2 c n . Likewise, for a trapped Fermi gas in the unitary limi t, th e virial expansion of contact provides an 
excellent explanation for the experimental measurement at T > 0.5Tp [l07| . This part will be reviewed in Sec. III. 

(iii) The virial expansion functions for the single-particle spectral function and dynamic structure factor have been 
determined (75l . |7(| , to the second order in fugacity. These results enable an important qualitative understanding of 
recent exp erimental measurements on momentum- resolved rf-spectroscopy (46|.|47| and two-photon Bragg spectroscopy 
|l07l . Il08|, for a trapped Fermi gas in the unitary limit at temperature down to the superfluid transition. The 
determination of the third virial expansion functions is straightforward, but involves much heavier numerical efforts. 
The virial expansion of dynamic structure factor and of single-particle spectral function will be reviewed in Sees. IV 
and V, respectively. 



II. VIRIAL EXPANSION OF EQUATION OF STATE 



Let us consider the virial expansion of the thermodynamic potential O, for a balanced spin-1/2 Fermi gas with 
equal spin populations (/i^ = /.t^ = fi). The spin-population imbalanced case with /i-j- fi± will be discussed at the 
end of the section. All the equations of state can be derived from the thermodynamic potential. By Taylor-expanding 
Q = —ksThiZ in the fugacity, where Z = 1 + zQi + z 2 Q 2 + ■ • • , the thermodynamic potential takes the form, 

O = -k B TQ 1 [z + b 2 z 2 + ■■■+ b n z n + •••], (17) 

where b n is referred to as the n-th (virial) expansion coefficient. Note that, by definition in Eq. ([5|) is given by 
ft( 3 ) = -k B TQ x b n . It is readily seen that, 

62 = {Qi - Qi/2) /Qi, (18) 
h = (Q 3 - Q1Q2 + Q?/3) /Qi, etc. (19) 

These equations give a general definition of virial expansion, which is applicable to both homogeneous and trapped 
systems. The calculation of the n-th virial coefficient requires the input of cluster partition function Qi (i < n), 
and hence requires the solutions of up to the n-particle problem. In practice, it is convenient to concentrate on 
the interaction effects only. We therefore consider the diffe rence A6„ = b n — bn and AQ n = Q n — Qn , where 
the superscript "1" denotes the non-interacting systems |l55| . For the second and third virial coefficients, one shall 
calculate respectively 

Ab 2 = AQ2/Q1 (20) 

and 



A6 3 = AQ3/Q1 - AQ 2 . (21) 

As we mentioned earlier, the calculation of virial coefficients in the strongly correlated regime is a subtle theoretical 
problem. The second virial coefficient was known long time ago through the elegant Beth-Uhlenbeck formalism, which 
relates in a s impl e manner the second virial coefficient to the two-body S'-matrix or the two-body scattering phase 
shift (97l. l98l Tl 1 1| . A connection between the virial series and the scattering matrix has been suspected since then. 
However, the computation of the third virial coe fficient, a long the line of Beth and Uhlenbeck's original work, met 
the very difficulties of the three-particle problem |l!2l . Ill3| . 

Until very recently, there is renewed interest in calculating higher-order virial coefficients, largely due to the cre- 
ation of ultracold atomic Fermi gases. Initial attempt was based on the field theoretic method, by calculating the 
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contribution of three- particle scattering proc ess t o the thermodynamic potential |ll4 |. It was shown by Rupak in 
2007 that in the unitary limit, Ab 3 ~ +1.05 [HI- 

However, it was soon realized by Liu, Hu and Drummond [72] that 
this value does not agree with the high-temperature heat-capacity measurements reported by Thomas's group at Duke 
University [2!| . Using an entire different strategy based on Eq. (|2"Tj) and the exact three-particle solution in harmonic 
traps, they predicted Ab% = —0.35510298. The numerical accuracy of A63 can be improved by including more three- 
particle energy levels. The latest calculation by Rakshit, Daily and Blume, along the line of Liu, Hu and Drummond's 
work, gave Ab 3 = -0.3551030264897 and A64 = -0.016 ± 0.004 In parallel, new field theoretic calculations for 
the third virial coefficient have been p erform ed. It was shown by Kaplan and Sun A63 = — 0.3573 ± 0.0005 |ll6j and 
by Leyronas A63 = —0.3551 ± 0.0001 |ll0| . At this stage, more complete field theoretic calculation is desirable, in 
order to confirm independently the fourth virial coefficient and to predict new virial coefficients. 

On the other hand, the experimental accuracy in measuring the equation of state of a unitary Fermi gas is improved 
very rapidly. The measurement by Salomon's group at Ecole Normale Superieure (ENS) gave Ab 3 ^ exp t = —0.35 ±0.02 
and Ap4 expt = 0.096 ± 0.015 [3(|. The latest measurement by Zwierlein's group at MIT reported Ab^ exp t = 0.096 ± 
0.010 J32J. We anticipate that the new predictions on the virial coefficients, improved continuously by many theorists, 
will play an important role in deepening our understanding of the equation of state of strongly correlated Fermi 
systems. 



A. Virial coefficients of non-interacting Fermi gases 



The background non-interacting virial coefficients can be conveniently determined by the non-interacting thermo- 
dynamic potential. For a homogeneous two-component Fermi gas, it takes the form [95|, 



-2k B T 



EM 1 



-(e k -ti)/(k B T) 



(22) 



where = h 2 k 2 /(2m) is the single-particle energy, the factor of 2 accounts for the spin degree of freedom. By using 
^ k = V f Q 47rfc 2 dfc/(27r) 3 and introducing a new variable t = e^/ (k B T), the ideal thermodynamic potential becomes, 



-V 



2k B T 2 



dt, 



(23) 



where \j B = [2irh 2 / (mksT)] 1 ' 2 is the thermal de Broglie wavelength. It is easy to identify Qi = 2V/\^ B . Therefore, 
by Taylor-expanding ln(l + ze _t ) in fugacity z and integrating out t term by term, we obtain the non-interacting 
virial coefficients in free space, 



b {1 



-1) 



71+1 



,5/2 



(24) 



For a Fermi gas in a harmonic trapping potential Vt(x) — muj^x 2 + y 2 + z 2 )/2, it is convenient to use the 
semiclassical approximation, or the so-called local density approximation. In the non- interacting limit, this amounts 
to setting, 



(i) 



r/x 



(x) 



v 



-2k B T 



-[ek+V T (x)-/i]/(fc B T) 



(25) 



where locally the single-particle energy is given by ek + Vr(x). Hereafter, we take the subscript "T" to denote the 
quantity in the trapped system, otherwise, by default we refer to a homogeneous system. As before, the integrations 
over x and k can be done by introducing a new variable t = [e^ + Vr(x)]/(fcsT) ■ This leads to, 



in 



2(k B Tf 

(IVjJt) 



i / t 2 In (l + ze~ l ) dt 



(26) 



where Qi,t = 2 (k B T) 3 / (Swr) 3 . By Taylor-expanding the log-term in fugacity z. we find the non-interacting virial 
coefficients in harmonic traps, 



b {1) 



-1) 



71+1 



(27) 
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We note that the use of semi-classical approximation means to neglect the discreteness of the energy spectrum in 
traps. Mathematically, this is equivalent to take a small parameter Cjt = fu^T/ikgT) and to keep in the results the 
leading term in lit- We note also that the non- interacting virial coefficients in the homogeneous case and trapped 
case are related by, b^ T = n~ a / 2 bn (n = 1, 2, 3, • • • ). 



B. Universal relation between homogeneous and trapped virial coefficients 

The correspondence relation discussed in Sec. IIA holds for a strongly interacting Fermi gas as well. Here, the 
crucial point is that the virial coefficients become temperature independent. To understand this, we note that in 
general the coefficients should be a function of the ratio XdB/a s , between the only two length scales XdB and a s - The 
temperature dependence enters through the thermal de Broglie wavelength. However, for a unitary Fermi gas where 
Ads/fls = 0, this dependence disappears. This is indeed a manifestation of fermionic universality, shared by many 
quantum systems with strong short-range interactions. 

In the thermodynamic limit, let us consider the thermodynamic potential of a harmonically trapped Fermi gas in 
the local density approximation, 

QT= _2_jkBrf_ [z + b 2 + b 3 + ... ]= [jjM., (28) 

where the trapped virial coefficients b n ,T are to be determined and fi(x) is the local thermodynamic potential, 

O(x) = -V^f- [z (x) + b 2 z 2 (x) + b 3 z 3 (x) + ■ ■ • ] . (29) 

X dB 

Here, the local fugacity z (x) = exp[/x(x) /(JzbT)] = zexp[— Vr(x)/(fcsT)] is given by the local chemical potential 
/x(x) = fj, — Vt(x). Because of the temperature- independent (constant) virial coefficients b n , the spatial integration 
can be done explicitly. This immediately leads to the universal relation for the virial coefficients of a unitary Fermi 

gas, 

<t = (30) 



C. Second virial coefficient of interacting Fermi gases 

1. Beth-Uhlenbeck formalism 

As shown by Beth and Uhlenbeck in 1937 [9?]], the second virial coefficient can be expressed in terms of the phase 
shifts of a two-body scattering problem. For a spin-1/2 Fermi gas, it takes the form, 



(21 + 1) 



dk—J-e 

dk 



(31) 



where the first summation is over all the two-body bound states (with the energy E B ) and 8i (k) is the phase shift of 
the Z-th partial wave. The second virial coefficient therefore can be determined for arbitrary interatomic interactions. 
For a pedagogical explanation of the elegant Beth-Ulenbeck formalism, we refer to the classical book by Kerson Huang 
[94[. This formalism has been applied by Ho and Mueller to explore the universal properties of atomic gases near 
a Feshbach resonance at high temperatures [70]. It has also been used extensively to study the equation of state of 
nuclear and neutron matter [87H92l |. 

For a s-wave Feshbach resonance that is of interest in ultracold atoms, the general expression for the s-wave phase 
shift Si =0 (k) is [H, & 



fccotJo (k) = h \r k 2 

a, 2 



which leads to, 



dSp 
dk 



r k 2 



r k' 



(32) 



(33) 
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Here for general discussion we have kept a nonzero range of interactions, ro- It is easy to see that the main contribution 
to the integral in Eq. (j3"Tj) , I = J °° dk(dSo/dk) exp(— \ 2 dB k 2 /2-k), comes from the region k ~ l/|a s |. Thus, after 
introducing a new variable y = k \a s \, the integral becomes, 



l + y 2 r /(2a s ) ( A 



/ = -sgn(a s ) / d» - '- y - ^ , exp ( -^y 2 ) . 



'[l-y 2 r /(2a s )] 2 + y 2 ~"" \ 2-no 2 ' 

U 

In the case of zero-range approximation (ro = 0), we obtain [70| . 

/(°) = -sgn(a s )|[l-erf(x)]e K2 , (35) 

where x = Xds/ (\/27r [ <x s | ) & n d erf (x) is the error function. The correction due to a finite range of interactions ro 
can be taken into account by Taylor-expanding the function in Eq. Q34p in a series of ro/ (2a s ). To the leading order, 
we find that, 



= ( ^ + 2^ 2 [1 - erf (*)] ) . 

4|o a | la; J 



(36) 



Near a Feshbach resonance where x -C 1 and ro -C |a s |, we have, 

In terms of the small dimensionless interaction strength {l/{kp a s ) <§; 1) and the range of interactions (kpro <C 1), 
the second virial coefficient can be written as, 

i— rp sgn(a,) 2 Tp 1 1 / T , 

V2 v 71 " » 2 fcfa s 4vtt V -if 

where the single bound state exists only for a positive scattering length with its energy i?e depending on both a s and 
ro. In the unitary limit, where l/(kpa s ) = 0, kpro = 0, and Eb = 0, we obtain the well-known result [94l. |98|. 

A0 2 - ^. (39) 

Concerning the experimental measurement, as an example, we estimate the second virial coefficient for 6 Li atoms 
using realistic experimental parameters. Let us consider the negative scattering length (BCS) side of the Feshbach 
resonance, for which the second virial coefficient takes the form, 

Ab 2 (a s <0) = ^ + ^=J^f-^—--^ x [^k F r . (40) 
V2 V n V 2 k F a s Ay/n V ±f 

The second and third terms on the right-hand-side of the above equation are non-universal since both of them depend 
on the temperature. These non-universal corrections are caused by a finite scattering length or a finite range of 
interactions. 

For fLi atoms, the finite scattering length near the Feshbach resonance Bq ~ 834 G can be conveniently calculated 
using (lli|, 

a s = a bg (l - ^L) , (41) 

where a^g ~ — 1405as in units of the Bohr radius as — 0.529 x 10 -10 m, and AB ~ 300 G. At the typical experimental 
density, where 1/kp ~ 400 nm, we find that kpa s ~ ±100, if the magnetic field is tuned away from the resonance 
by one Gauss. This leads to about a percent correction to the second virial coefficient at the d egen erate temperature 
Tp. On the other hand, the finite range of interactions near the resonance can be modeled as |119j . 

ro= - 2a ( 1 _f%y + «-^, ,42, 

where i?„ ~ 0.0269 nm and b ~ 2.1 nm is essentially the Van der Waals length. As i?* <C b <C |a s | across the Feshbach 
resonance, the finite range of interactions is reduced to a constant ro ~ 4.7 nm. Thus, we obtain the dimensionless 
range of interactions kpr a ~ 0.012, for the typical Fermi wavelength. It gives about 0.1% correction to the second 
virial coefficient at Tp. 
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2. Field theoretic method 



The second virial coefficient can also be conveniently calculated using the field theoretic method |117l |. i.e., the 
diagrammatic expansion method we mentioned earlier in Sec.lB. This provide a simple example to illustrate the close 
relation between the virial expansion and the diagrammatic expansion. In the following, we introduce briefly the 
procedure. To obtain the virial coefficients, the basic idea of field theoretic meth od is to calculate &S l \ which is the 
contribution of /-particle scattering process to the thermodynamic potential [l!5| . At large temperatures, we expand 
in fugacity, 

^>=-F^f>«z", (43) 

where b$ is the n-th virial coefficient from /-particle interaction. The total n-th virial coefficient b n = b n + b^ + 
■•• + bn , where b n ^ = (— l) n+1 n~ 5 / 2 is the virial coefficient of an ideal Fermi gas. In the case of I = 2, we have 

Ab 2 = & 2 -4 1) = 4 2) - 

We start from a path-integral functional action (l9l . [53| , using the single-channel fermionic model with zero-range 
interactions Eq. (|10[) . By performing a Hubbard-Stratonovich transformation to decouple the interaction term, the 
original fermionic partition function Z — J T>[ip(x), ■0(x)]e^ s can be expressed as Z = J D[A(x), A* (x)]e~ s "ff , in 

terms of bosonic variables A(x). The "effective" bosonic action can be written in a series expansion: S e ff = ^eff 
In the normal state, the first term in the expansion reads [l||, Ql| [53[ , 

^^[-xWlA^A'lg), (44) 



where 



4irK 2 a s V ^ 

k 



/iK£q/2+k) + /F(Cq/2-k) ~ 1 \_ 

iv n - £q/2+k - Cq/2-k 2 £ k 



(45) 



is the two-particle propagator. Here we have used the abbreviation q = (c\,iv n ), the bosonic (fermionic) Matsubara 
frequency v n = 2nirkBT (oj m = (2m + 1)7t/sbT), £k = £k — " = fi 2 k 2 /(2m) — u, and the Fermi distribution function 
Jf (x) = l/[exp(x/fcsT) + 1]. The action S e ^ accounts for the scatterings between two-particles in the presence of 

('2) 

other particles (i.e., medium), and thus includes the two-body contribution to all the virial coefficients of b„ (n > 2). 
It gives rise to the following thermodynamic potential given by Nozieres and Schmitt-Rink (NSR) in 1985 (la ]. 



1 f +oc 

^ = kBT J2 In [-x (q)} cxp (iv n 0+) = -~J2 dQ fs (n) 5 (q, fi) , (46) 

where the summation over the Matsubara frequency has been converted into an integral using a phase shift, 

8 (q, 0) = - Imln [- X (q, iv n -> fi+)] , (47) 

and fs (x) = l/[exp(cc/fcsT) — 1] is the Bose-Einstein distribution function. The diagrammatic representation of fl^ 
has been illustrated earlier in Fig. 1. 

In the high-temperature limit, where the fugacity z = exp^/kgT) <C 1 and ./f (£k) ~ zexp(— ^k/^s? 1 ), we 
may Taylor-expand the phase shift in powers of z. Focusing on the unitary limit, we approximate the two-particle 
propagator 

X (q, 0+) = X (0) (q, 0+) + z X {1) (q, 0+) + 0(z 2 ), (48) 

where 

(i) \ - exp (-C q /2+k/fci?r) + exp (-jq/g-k/fcgT) , s 

X ^ fi+ - e q /2 - 2e k + 2u ' 1 J 
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and fl + = O + i0 + . To the leading order of x' ' (q, the phase shift is exactly a step function, 

5W(q,0) = |e(f2-^+2 M ). (51) 
Thus, to the leading order of fugacity we have, 

This gives rise to the second virial coefficient b% = A&2 = 1/ V2- Away from the unitary limit, it is straightforward to 
show that we can recover the Beth-Uhlenbeck formalism from Eq. (|46p . by taking the phase shift (5(q, S7) in vacuum. 

We note that the higher-order contribution of fe„ ; (n > 3) can be obtained by successively calculating the z n term 
inEq. (gU). 



D. Virial coefficients from exact few-body solutions in harmonic traps 

We now turn to calculate the third virial coefficient, by using an entirely different method [73}. We solve first 
the two-particle and three-particle problems in an isotropic 3D harmonic trap Vr(x) = muj 2 ^x 2 /2, and then use the 
solutions to obtain the second and third virial coefficients. In the end, we discuss the possibility of calculating the 
fourth virial coefficient. 

We note that in cold-atom experiments the harmonic trap is often highly anisotropic. The three-particle problem 
at unitaritiy in an anisotropic trap can hardly be solved exactly. Fortunately, for a large number of particles, for 
which the local density approximation is valid, we are free to use harmonic traps of any aspect ratio to calculate the 
virial coefficients, by using the universal relation Eq. (|50|) . 



1. Relative Hamiltonian of few-particle systems 



In a harmonic trap, it is useful to separate the center-of-mass motion and relative motion. We thus defin e th e 
following center-of-mass coordinate R and relative coordinates (i > 2) for N fermions in a harmonic trap (99l . llOfjj . 

R= (xi + ---+Xiv)/JV, (53) 

and 

r ' = /¥( x -rig4 (54) 

respectively. In this Jacobi coordinate, the Hamiltonian of the non-interacting Schrodinger equation takes the form 
Ho = H cm + Ureh where, 

n cm = -^l+lM^R 2 , (55) 

and 



N 

E 



i 



+ -m<4r, 2 



(56) 



The center-of-mass motion is simply that of a harmonically trapped particle of mass M = Nm, with well-known wave 
functions and spectrum E cm = (n cm + 3/2)huTi where n cm = 0, 1,2... is a non-negative integer. In the presence of 
interactions, the relative Hamiltonian should be solved in conjunction with the Bethe-Peierls boundary condition, Eq. 
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2. Two fermions in a 3D harmonic trap 
Let us consider the two-fermion problem in a harmonic trap, where the relative Schrodinger equation becomes 

h 2 _, 1 



^t{v)=E rel ^t(v), (57) 



where two fermions with unlike spins do not stay at the same position (r > 0). Here, we have re-defined r = \/2r2 and 
without confusing with the chemical potential we have used a reduced mass /i = m/2. It is clear that only the I = 
subspace of the relative wave function is affected by the s-wave contact interaction. According to the Bethe-Peierls 
boundary condition, as r — > the relative wave function should take the form, ip^fir) (l/ r — l/ a s)j or satisfy, 
d[rtp2b j l® r = — (j"^2b l ) / a s- The two-fermion problem in a harmonic trap was first solved by Busch and coworkers 
[80(. In the following, we present a simple physical interpretation of the solution. 

The key point is that, regardless of the boundary condition, there are two types of general solutions of the relative 
Schrodinger equation ([ST)) in the I = subspace, ( r ) ^ exp(— r 2 /2d 2 )/(r/<i). Here the function f(x) can either 
be the first kind of Kummer confluent hypergeometric function %Fi or the second kind of Kummer confluent hyper- 
geometric function U. We have taken d = \J h/JjuJr) as the characteristic length scale of the trap. In the absence 
of interactions, the first Kummer function gives rise to the standard wave function of 3D harmonic oscillators. With 
interactions, however, we have to choose the second Kummer function U, since it diverges as 1/r at origin and thus 
satisfies the Bethe-Peierls boundary condition. 

Therefore, the (un-normalized) relative wave function and relative energy should be rewritten as, 

<(r;^) = r(-i/)tf(-i/,|,^)exp(-^), (58) 

and 

E re i = {2v + |)fiwr, (59) 

respectively. Here, T is the Gamma function, the real number v plays the role of a quantum number and should be 
determined by the boundary condition, lim r _>o d (?"02&') /dr = — (r^f) jo. By examining the short range behavior 
of the second Kummer function U(— u, 3/2, x), this leads to the familiar equation for energy levels [8C| . 



2T(-v) d 



r(-i/-l/2) a s 



(60) 



In Fig. [5J we give the resulting energy spectrum as a function of the dimensionless interaction strength d/a s . 

The spectrum is easy to understand. At infinitely small scattering length a s — > CP, v(a s = 0~) = n re i (n re i = 
0,1,2...), which recovers the spectrum in the non-interacting limit. With increasingly attractive interactions, the 
energies decrease. In the unitarity (resonance) limit where the scattering length diverges, a s — > ±oo, we find that 
v(a = ±oo ) = n re i — 1/2. As the attraction increases further, the scattering length becomes positive and decreases 
in magnitude. We then observe two distinct types of behavior: the ground state is a molecule of size a, whose energy 
diverges asymptotically as —h 2 /ma 2 s as a s — > + , while the excited states may be viewed as two repulsively interacting 
fermions with the same scattering length a s . Their energies decrease to the non-interacting values as a s — > + . 

In this two-body picture, a universal repulsively interacting Fermi gas with zero-range interaction potentials may be 
realized on the positive scattering length side of a Feshbach resonance for an attractive interaction potential, provided 
that all two fermions with unlike spins occupy the exited states or the upper branch of the two-body energy spectrum. 



3. Three fermions in a 3D harmonic trap: General exact solutions 



Let us turn to the three fermion case by considering two spin- up fermions an d on e spin-down fermion, i.e., the 
configuration shown in Fig. |H1 The relative Hamiltonian can be written as (99l . Il00| , 

Urel - ^ (V r 2 + Vl) + (r 2 + p 2 ) , (61) 

where we have redefined the Jacobi coordinates r = s/2r2 and p = ^/2v 3 . which measure the distance between the 
particle 1 and 2 (i.e., pair), and the distance from the particle 3 to the center-of-mass of the pair, respectively. 
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Figure 5: (color online) Energy spectrum of the relative motion of a trapped two-fermion system near a Feshbach resonance 
(i.e, d/a s — 0, where d is the characteristic harmonic oscillator length). For a positive scattering length a s > in the right part 
of the figure, the ground state is a molecule with size a B , whose energy diverges as E re i ~ — ft 2 / (m«J). The excited states or the 
upper branch of the resonance may be viewed as the Hilbert space of a "repulsive" Fermi gas with the same scattering length a s . 
In this two-body picture, the level from the point 2 to 3 is the ground state energy level of the repulsive two-fermion sub-space, 
whose energy initially increases linearly with increasing a a from 1.5/kjT at the point 2 and finally saturates towards 2.5/kJT at 
the resonance point 3. For comparison, we illustrate as well the ground state energy level in the case of a negative scattering 
length and show how the energy increases with increased scattering length from point 1 to 2. From ref. (73| ; copyright (2010) 
by APS. 




Figure 6: (color online) Configuration of three interacting fermions, two spin-up and one spin-down. From ref. [73l |: copyright 
(2010) by APS. 

Inspired by the two-fermion solution, it is readily seen that the relative wave function of the Hamiltonian (|6ip may 
be expanded into products of two Kummer confluent hypergeometric functions. Intuitively, we may write down the 
following ansatz |72| . 

^(r,p) = (l-P 13 )x(r,p), (62) 

where, 

x (r, P ) = J2 ^ r 2b(r; n, n )Rni (p) rr (p) ■ (63) 

The two-body relative wave function r (p2b l { r ^ v i,n) with energy (2i/;„ + 3/2)hu>T describes the motion of the paired 
particles 1 and 2, and the wave function R n i (p) Y™ (p) with energy (2n + / + 3/2)hidT gives the motion of particle 
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3 relative to the pair. Here, i? n j (p) is the standard radial wave function of a 3D harmonic oscillator and Yj m (p) is 
the spherical harmonic. Owing to the rotational symmetry of the relative Hamiltonian (|6ip . it is easy to see that the 
relative angular momenta I and m are good quantum numbers. The value of i>i m is uniquely determined from energy 
conservation, 



Erei = [(2«j, n + 3/2) + (2n + I + 3/2)] fiw T , 



(64) 



for a given relative energy E re i. It varies with the index n at a given angular momentum I. Finally, V\3 is an exchange 
operator for particles 1 and 3, which ensures the correct exchange symmetry of the relative wave function due to Fermi 
exclusion principle, i.e., Vv&x{ r i p) = X ( r /2 + V^p/2, V3r/2 — p/2). The relative energy E re i together with the 
expansion coefficient a n should be determined by the Bethe-Peierls boundary condition, i.e., lim r _>.o [drip^ f (r, p)]/dr = 
— [riplf (r, p)]/a s . We note that the second Bethe-Peierls boundary condition in case of particle 2 approaching particle 
3 is satisfied automatically due to the exchange operator acting on the relative wave function. 
By writing \ ( r ; P) = <K r J PW™ {p)i the Bethe-Peierls boundary condition takes the form (r->0), 



— [r(f>(r, p}] 
a. 



d[r<f>(r,p)] 
dr 



(-1)' 



. Vip p^ 

■ 2 ' 2' 



Using the asymptotic behavior of the second kind of Kummer function, lim^ 
2y/TrT (—Vi n) /r (—v^n — 1/2), it is easy to show that in the limit of r — > 0, 



,.0 r {-Vl, n ) U(-Vl. 



[r<f>(r, p)} 



/7F \ ^ 

— > a n R n i (p) . 
1 - ^ — ^ 



(65) 

, 3/2, x 2 ) =^/x- 
(66) 



and 



d [r<p(r, p)} 
dr 



= -\/Tr^a n R n i (p) 



2r(-n,„) 
r(-i^, n -i/2)' 



Thus, the Bethe-Peierls boundary condition becomes, 



B n R nl (p)-R nl [- 



where 



d 
a s 



2T{-v hn ) 



= 0, 



r(-^,„-l/2) 

Projecting onto the orthogonal and complete set of basis functions R n i (p), we find that a secular equation, 

ZT(-vi, n ) ^ „ f d 

T{-v Ln -l/2) an + 

where we have defined the matrix coefficient, 



00 

C nn > = J p 2 dpR rd (p) R n > 



(67) 



(68) 



(69) 



(70) 



(71) 



which arises from the exchange effect due to the operator V13. In the absence of C nn i , the above secular equation 
describes a three-fermion problem of a pair and a single particle, un-correlated to each other. It then simply reduces 
to Eq. (|60|) . as expected. 

The secular equatio n (1701 was first obtained by Kestner and Duan by solving the three-particle scattering problem 
using Green function J82j]. To solve it, for a given scattering length we may try different values of relative energy 
E re i, implicit via v^ n - However, it turns out to be more convenient to diagonalize the matrix A = {A„„/} for a given 
relative energy, where 



A , - 



2r(-n, n ) 

1/2) n "' 



(72) 



The eigenvalues of the matrix A then gives all the possible values of d/a s for a particular relative energy. We finally 
invert a(E re i) to obtain the relative energy as a function of the scattering length. Numerically, we find that the matrix 
A is symmetric and thus the standard diagonalization algorithm can be used. We outline the details of the numerical 
calculation of Eq. (|T2"|) in the Appendix A. 
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4- Three fermions in a 3D harmonic trap: Exact solutions in the unitarity limit 



In the unitarity limit with infinitely large scattering length , a. — > o o, we may obtain more physical solutions using 
hyperspherical coordinates, as shown by Werner and Castin |8ll . Il0dl |. By defining a hyperradius R = y 1 (r 2 + p 2 )/2 

and hyper angles fl = (a, f , p), where a — arctan(r/p) and r and p are respectively the unit vector along r and p, we 
may write (suHooj, 



)=^(1-^3)^^(P), (73) 
/ R sm (2a) 



to decouple the motion in the hyperradius and hyp erangles for given relative angular momenta I and m . It leads to 
the following decoupled Schrodinger equations |l00l |. 



1 , $ 



and 



F" -ji F '+[ j£+ <4^ 2 F = 2E relF, (74) 



v " (a)+ l l^±y (a) = (a), (75) 



where Sj is the eigenvalue for the n-th wave function of the hyperangle equation. 

For three- fermions, sf n is always positive. Therefore, the hyperradius equation (|74p can be interpreted as a 
Schrodinger equation for a fictitious particle of mass unity moving in two dimens ions in an effective potential 
(sf n /R 2 + Ld^R 2 ) with a bounded wave function F(R). The resulting spectrum is [HI, llOOf 

E rel = (2q + si <n + 1) hw T , (76) 

where the good quantum number q labels the number of nodes in the hyperradius wave function. 

The eigenvalue s;.„ shou ld be determined by the Bethe-Peierls boundary condition, which in hyperspherical coor- 
dinates takes the from (sT|, Il00| , 

tf(0)-(-iy± <p (l)=Q. (77) 



In addition, we need to impose the boundary condition tp (tt/2) = 0, since the relative wave function (|73[) should not 
be singular at a = tt/2. The general solution of the hyperangle equation (|75[) satisfying ip (tt/2) = is given by, 

-1 IP fl + l-Sl,n l+l + SLn , 3 _ 2 



<P « x^ 2 F, ^ ~ 2 ' 2 + r ; ^ j , (78) 

where a; = cos(a) and 2-F1 is the hypergeometric function. In the absence of interactions, the Bethe-Peierls boundary 
condition ([77]) should be replaced by ip (0) = 0, since the relative wave function ([73)) should not be singular at 
a = either. As <p(Q) = T(l + 3/2)T(l/2)/[r((Z + 2 + s i>n )/2)r((Z + 2 - s;,„)/2)], this boundary condition leads 

to [I + 2 — s[„]/2 = — n, or sj 1 ^ = 2n + Z + 2, where n. = 0, 1, 2, ... is a non-negative integer and we have used the 
superscript "1" to denote a non- interacting system. However, a spurious solution occurs when I = and n = 0, for 
which sj 1 ^ = 2, v?(a) = sin(2a)/2 and thus, the symmetry operator (1 — V13) gives a vanishing relative wave function 



in Eq. (|73|) that should be discarded [100]. Wo conclude that for three non-interacting fermions, 

(i) f2n + 4, Z = 

~ \ 2n + l + 2, I > ' [ ' 

For three interacting fermions, we need to determine s^ n by substituting the general solution (|78[) into the Bethe- 
Peierls boundary condition (|77l) . In the Appendix B, we describe how to accurately calculate s; in . In the boundary 
condition Eq. (]77| . the leading effect of interactions is carried by tp' (0) and therefore, tp' (0) = determines the 
asymptotic values of s; iM at large momentum Z or n. This gives rise to (Z + 1 — s^ n )/2 = —n, or, 

/ 2n + 3, Z = , an , 

S ''"- \ 2n + Z + l, Z>0 ' (80) 

where we have used a bar to indicate the asymptotic results. By comparing Eqs. ([75]) and (]5Dj) . asymptotically the 
attractive interaction will reduce s/ jn by a unity. 
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5. Three fermions in a 3D harmonic trap: Energy spectrum 



We can numerically solve both the general exact solution (f62j) along the BEC-BCS crossover and the exact solution 
(|73p in the unitarity limit. In the latter unitary case, the accuracy of results can be improved to arbitrary precision 
by using suitable mathematical software, described in Appendix B. Fig. [7] reports the energy spectrum of three 
interacting fermions with increasingly attractive interaction strength at the ground state angular momentum, 1 = 1. 
For a given scattering length, we typically calculate several ten thousand energy levels (i.e., E re i < (I + 256)/kjy) 
in different subspace. To construct the matrix A, Eq. (|72j) . we have kept a maximum value of n max = 128 in the 
functions R n i (p). Using the accurate spectrum in the unitarity limit as a benchmark, we estimate that the typical 
relative numerical error of energy levels is less than 10~ 6 . We have found a number of nontrivial features in the energy 
spectrum. 



Figure 7: (color online) Relative energy spectrum of three interacting fermions at the ground state subspace 1 = 1. On the 
positive scattering length (BEC) side of the resonance, there are two types of energy levels: one (is vertical and) diverges with 
decreasing the scattering length a s and the other (is horizontal) converges to the non-interacting spectrum. The latter may be 
viewed as the energy spectrum of three repulsively interacting fermions. In analogy with the two-fermion case, we show the 
ground state energy level of the repulsive three- fermion system (i.e. from point 2 to 3), as well as the ground state energy level 
of the attractive three-fermion system for a s < (i.e., from the point 1 to 2). In the unitarity limit, we show by the circles the 
energy levels that should be excluded when we identify the energy spectrum for infinitely large repulsive interactions. Adapted 
from ref. [zl; copyright (2010) by APS. 

The spectrum on the BCS side is relatively simple. It can be understood as a non-interacting spectrum at d/a s —> 
— oo, in which E re i = (2Q + 3)Hujt at Z = and E re i = (2Q + 1 + 1)Hujt at I > 1, with a positive integer Q = 1, 2, 3, ... 
that denotes also the degeneracy of the energy levels. The attractive interactions reduce the energies and at the same 
time lift the degeneracy. Above the resonance or unitary point of d/a s = 0, however, the spectrum becomes much 
more complicated. 

There are a group of nearly vertical energy levels that diverge towards the BEC limit of d/a s — > +oo. From the 
two-body relative energy spectrum in Fig. [5l we may identify these as energy states containing a molecule of size a s 
and a fermion. For a given scattering length, these nearly vertical energy level differ by about 2/iwt, resulting from 
the motion of the fermion relative to the molecule. In addition to the nearly vertical energy levels, most interestingly, 
we observe also some nearly horizontal energy levels, which converge to the non-interacting spectrum in the BEC 
limit. In analogy with the two-body case, we may identify these horizontal levels as the energy spectrum of three 
repulsively interacting fermions. We show explicitly in the figure the ground state level of three repulsively interacting 
fermions, which increases in energy from the point 2 to 3 with increasing scattering length from a s = + to a s = +oo. 
For comparison, we also show the ground state level of three attractively interacting fermions at a negative scattering 
length, which decreases in energy from the point 2 to 1 with increasing absolute value of a s . 

This identification of energy spectrum for repulsive interactions, however, is not as rigorous as in the two-body case. 
There are many apparent avoided crossings between the vertical and horizontal energy levels. Therefore, by changing a 
positive scattering length from the BEC limit to the unitarity limit, three fermions initially at the horizontal level may 
finally transition into a vertical level, provided that the sweep of scattering length is sufficiently slow and adiabatic. 
This leads to the conversion of fermionic atoms to bosonic molecules. A detailed analysis of the loss rate of fermionic 




10 



20 



atoms as a function of sweep rate may be straightforward obtained by applying the Landau-Zener tunnelling model. 

Let us now focus on the resonance case of most significant interest. In Fig. [JJ we show explicitly by green dots 
the vertical energy levels in the unitarity limit. These levels should be excluded if we are interested in the spectrum 
of repulsively interacting fermions. Amazin gly, for each given angular momentum, these energy levels form a regular 
ladder with an exact energy spacing 2hujT [99(- Using the exact solution in the unitarity limit, Eq. (|76p . we may 
identify unambiguously that the energy ladder is given by, 

E ret = (2q + si t0 + l)hujT. (81) 

Therefore, in the unitarity limit the lowest-order solution of the hyperangle equation gives rise to the relative wave 
function of a molecule and a fermion. Thus, it should be discarded when considering three resonantly interacting 
fermions with an effective repulsive interaction. 



6. Second virial coefficient 



We now calculate the virial coefficients of a trapped attractively interacting Fermi gas. In a harmonic trap, the 
oscillator length d provides a large length scale, compared to the thermal wavelength XdB ■ Alternatively, we may use 
lot = fi^T / '{ksT) < 1 to characterize the intrinsic length scale relative to the trap. All the virial coefficients and 
cluster partition functions in harmonic traps therefore depend on the small parameter &t- We shall be interested in 
a universal regime with vanishing Cup, in accord with the large number of atoms in a real experiment. 

To obtain A&2,Tj we consider separately AQ2,t and Qi,t- The single-particle partition function Qi,t is de- 
termined by the single-particle spectrum of a 3D harmonic oscillator, E n \ = (2n + I + 3/2)fiwr- We find that 
Qi,t = 2/[exp(+tJT/2) — exp(— wt/2)] 3 — 2 (ksT) / (hiop) , in agreement with the previous result based on the local 
density approximation (see Eq. |21))) . The pre- factor of two accounts for the two possible spin states of a single fermion. 
In the calculation of AQ2.T, it is easy to see that the summation over the center-of-mass energy gives exactly Qi.p/2. 
Using Eq. ([59)1 . we find that, 



Ab 



2,T 



where the non-interacting Vn = n (n — 0, 1, 2, ...). 

At resonance with an infinitely large scattering length, the spectrum is known exactly: v n = n — 1/2, giving rise to, 

r c^h-^^ 1_ 1_ 



_ 1 cxp(-o; T /2) 
2,T ~ 2 [1 + exp (-U) T )} 



The term ujj, in Eq. (|83p is nonuniversal and is negligibly small for a cloud with a large number of atoms. We 
therefore obtain the universal second virial coefficient: A&2,t = 1/4, which are temperature independent. 

In Fig. [5) we show the second virial coefficient through the BEC-BCS crossover at three typical temperatures. Here 
we consider a gas with N = 100 atoms and scale the inverse scattering length using the Fermi vector at the trap center, 
kp = (24A r ) 1 / 6 /(d/v / 2)- The temperature is given in units of Fermi temperature Tp = Ep/ks = (3N) 1 ^ 3 (hujT /ks)- 
All the curves with distinct temperatures cross at a s — > ±00. This is the manifestation of universal behavior anticipated 
if there is no any intrinsic length scale. However, the characteristic length scale d of harmonic traps brings a small 
(non- universal) temperature dependence that decreases as N~ 2 / 3 , shown by the terms £j\ in Eq. ([8"3")h 

According to the universal relation between trapped and homogeneous virial coefficients, Eq. (|30[) . we obtain 
immediately the homogeneous second virial coefficient in the unitarity limit, A62 = l/y/2, which is in agreement with 
the result obtained from the Beth-Uhlenbeck formalism and from the field theoretic calculation. 



7. Third virial coefficient 



The calculation of the third virial coefficient, which is given by Ab 3 t = &-Qs,t /Qit — AQ 2 ,t, is more complicated. 
Either the term AQ?,^ / Qi.t or AQ2,t diverges as ojt —> 0, but the leading divergences cancel with each other. In 
the numerical calculation, we have to carefully separate the leading divergent term and calculate them analytically. 
It is readily seen that the spin states of an d configurations contribute equally to Q% : t- The term Qi t T in the 
denominators is canceled exactly by the summation over the center-of-mass energy. We thus have 



AQ 3) r/Qi,T = Y, c M-Erei/k B T) - £ cxp(-^ J ) /fc B T) . 



(84) 
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Figure 8: (color online) Second virial coefficient of a trapped attractive Fermi gas as a function of the interaction parameter 
l/(fcj-a s ). We have used a total number of atoms N = 100, leading to Cjt = (3-/V)~ 1//3 « 0.15 at T = Tf- Adapted from ref. 
El; copyright (2009) by APS. 



To proceed, it is important to analyze analytically the behavior of E re i at high energies. For this purpose, we 
introduce a relative energy E re i , which is the solution of Eq. (|72p in the absence of the exchange term C' nm , and can 
be constructed directly from the two-body relative energy. In the subspace with a total relative momentum I, it takes 
the form, 



E re i = {2n + l + 3/2) huj T + {2v + 3/2)fiw T , 



(85) 



where v is the solution of the two-body spectrum of Eq. At high energies the full spectrum E re i approaches 

asymptotically to E re i as the exchange effect becomes increasingly insignificant. There is an important exception, 
however, occurring at zero total relative momentum 1 = 0. As mentioned earlier, the solution of E re i at n = and 
I = is spurious and does not match any solution of E re i. Therefore, for the I = subspace, we require n > 1 in Eq. 

It is easy to see that if we keep the spurious solution in the I = subspace, the difference [^2 cxp(— E re i/ksT) — 

^ exp(— E^/kgT)] is exactly equal to AQ 2 .t, since in Eq. (|85p the first part of spectrum is exactly identical to the 
spectrum of center-of-mass motion. The spurious solution gives a contribution, 



-(2v< 1) +3)iJT 



= 2e" 3 " T/2 A6 



2,T, 



(86) 



which should be subtracted. Keeping this in mind, we finally arrive at the following expression for the third virial 
coefficient of a trapped Fermi gas with attractive interactions: 



Ab, 



3,T 



E 



B — £ 



2e~ 3 ° T/2 Ab 



2.T- 



(87) 



The summation is over all the possible relative energy levels E re i and their asymptotic values E re i . It is well-behaved 
and converges at any scattering length. The third virial coefficient of a trapped attractive Fermi gas in the BEC-BCS 
crossover was shown in Fig. [9l 

In the unitarity limit, it is more convenient to use the exact spectrum given by Eq. (|76p . where s; >n can be obtained 

numerically to arbitrary accuracy and the non-interacting is given by Eq. ([79)1 . To control the divergence problem, 
we shall use the same strategy as before and to approach s;.„ by using its asymptotic value s^ n given in Eq. (|80[) . 
Integrating out the q degree of freedom and using Eq. ([83)) to calculate AQ2,t, we find that, 



Ab, 



3,T 



3 — 2u)T 



E 

/ . 1 1 



(88) 
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where A is given by 



„(!) 



A = > e~ UT ^" - e.- UTS '^ ? . (89) 

{.1 



(1-e 



We note that for the summation, implicitly there is a pre-factor (2/ + 1), accounting for the degeneracy of each 
subspace. The value of A can then be calculated analytically, leading to, 

A = -e-** (1 - e-**) . (90) 



We have calculated numerically J^i „( e " TSi ' n — e WtS! >») by imposing the cut-offs of n < n max = 512 and I < I 
512. We find that, 



max 



A6 3 , T ^ -0.06833960 + 0.038867w£ + • • • . (91) 

The numerical accuracy can be further improved by suitably enlarging n max and Z max . By neglecting the dependence 
on Q in the thermodynamic limit, we obtain the universal third virial coefficient: Ab 3t T — —0.06833960. Using 
the universal relation between trapped and homogeneous virial coefficients, Eq. (|30[) . we obtain immediately the 
homogeneous third virial coefficient in the unitarity limit, A63 ~ —0.35510298. 

In a recent study by Rakshit, Daily and Blume [78j], much more energy levels are included in the calculation of 
the third virial coefficient in the unitary limit. As a result, the accuracy is much improved. It was shown that [78j 
Ab 3T = -0.068339609311287 and Ab 3 = -0.3551030264897. 



8. Fourth virial coefficient 



The calculation of the fourth virial coefficient could follow the same strategy. However, the determination of AQ4 7- 
appears to be a daunting task, since so far the problem of four interacting fermions in harmonic traps has no exact 
solutions. 

This difficulty was overcome by Rakshit, Daily and Blume [78]], by using a scheme that allows to extrapolate the 
high temperature behavior of the virial coefficients from the low-lying portion of the excitation spectra only. This 
scheme is largely due to the weak Cjt dependence of the trapped virial coefficient A6„ t: because of the peculiarity of 
the harmonic trapping potential, Ab n ^T is a function of u)^. As a result, one can determine A6„.t at relatively large Cjt 
(i.e., lot ~ 1) and then extrapolate it to the zero-cJr limit. This procedure requires a small portion of the excitation 
spectra, which can be calculated using the stochastic variational approach (85| . with moderate computational resources. 
It was predicted that A6 4 ,t = -0.0020 ± 0.0005 and A64 = -0.016 ± 0.004. 

By using the same token, Rakshit, Daily and Blume estimated the fifth virial coefficient, 0.0017 < A65 < 0.101, 
and conjectured the sign of the higher-order virial coefficients is +, — , — , +,+,—,••• for n = 6, 7, 8, 9, 10, 11, • • • . 
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E. Third virial coefficient from field theoretic method 



Here we review briefly the diagrammatic calculation of the third virial coefficient. The basic idea is to calculate 
• t?( 3 ) = — dfl^ /dfi, wh ich involv es the contribution from the three-particle scattering process. As the three-particle 



or n 



vertex function is solved |f 2dl If 2 j |. in principle the third virial coefficient could be determined. However, as we shall 
see, the calculation turns out to be subtle. The diagrammatic representation of fl^ is shown in Fig. [TUJ The two- 
and three-particle vertex functions are indicated by T 2 and X3, respectively. 




Figure 10: (color online) Diagrammatic representation of the contribution of three-particle scattering process to the t herm ody- 
nam ic potential. Here T2 and T3 are respectively the two- and three-particle vertex functions. For details, see refs. |l2ll | and 

ma. 



The calculation of Q ^ at large temperatures was pioneered by Rupak |ll5j . by using a two-channel model for the 
description of Feshbach resonances. A dimer field is introduced, designed to reproduce the continuum two-body phase 
shift. In the unitary limit, it was predicted that A63 ~ 1.05. This calculation was recently improved by Kaplan and 
Sun [116], with the development of a new diagrammatic method for —cKl (A) /()(i. The sum over discrete Matsubara 
frequencies is converted to a Possion resummation. This leads to A 63 = — 0.3573 ± 0.0005. 
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Figure 11: (color online) Third virial coeffi cient as a function of the dimensionless parameter \dB/a„. The inset shows the 
second virial coefficient. Adapted from ref. [ll0l |. 



The latest field theoretical calculation of the third virial coefficient was given by Leyronas |ll0| , by using the single- 
channel Hamiltonian and Feynman diagrams for —dQ^/dfi or the single-particle Green function. Explicit analytic 
expressions of virial coefficient were obtained. To the accuracy of four digits, it was found that A&3 = — 0.3551±0.0001, 
which is in excellent agreement with the calculation based on the e xact three-particle solutions in harmonic traps 
but disagrees slightly with that obtained by Kaplan and Sun |ll6| . Fig. [Tl] shows the prediction by Leyronas 



11C 1 

At this stage, we believe that the result of A63 = — 0.3551±0.0001 is robust, as it has been checked independently by 
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two entirely different methods. The discrepancy between different field theoretic calculations remains to be unders tood . 
We note that it is appealing to calculate the fourth viri al co efficient A&4, along the line of Leyronas's calculation |llOl |. 
as the four-particle vertex function is basically known |l2l| . Together with an improved calculation with four-fermion 
solutions in harmonic traps, A64 could be determined very accurately. 



F. Virial equation of state for ultracold Fermi atoms and its comparison with experimental measurements 

1. Virial equation of state 

We are now ready to calculate the virial equations of states in the high temperature regime, by using the thermo- 
dynamic potential 

Q = fi« - V ^§^(Ab 2 z 2 + Ab 3 z 3 + ---) (92) 

*dB 

and 

Q T = Q« _ 2 ^ bT ^ (Ab 2 . T z 2 + Ab 3 . T z 3 + ■■■), (93) 
(Hlot) 

respectively, for a homogeneous or a harmonically trapped Fermi gas. Here, the non-interacting thermodynamic 
potentials are given by Eqs. (f2"3"| and (f2T))) . All the other thermodynamic quantities can be derived from the ther- 
modynamic potential by the standard thermodynamic relations, for example, N = — 9il/<9/i, S = —dQ/dT, and then 
E = n + TS + fxN, 

As an concrete example, let us focus on the unitary limit in the thermodynamic limit, which is of the greatest 
interest. The equations of states are easy to calculate because of the temperature independence of virial coefficients. 
It is also easy to check the well-known scaling relation in the unitarity limit: E = —30/2 for a homogeneous Fermi 
gas and — — 30 for a harmonically trapped Fermi gas [Sf3]. The difference of the factor of two arises from the 
fact (virial theorem) that in harmonic traps the internal energy is exactly equal to the trapping potential energy. 

To be dimensionless, we take the Fermi temperature Tp or Fermi energy (Ep = kpjTp) as the units for tem- 
perature and energy. For a homogeneous or a harmonically trapped Fermi gas, the Fermi energy is given by 
Ep = h 2 (3ir 2 N/V) 2 / 3 /2m and Ep = (SNy^hajp, respectively. In the actual calculations, we determine the number 
of atoms N, the total entropy 5, and the total energy E at given fugacity and a fixed temperature, and consequently 
obtain the Fermi temperature Tp and Fermi energy Ep. We then plot the energy or energy per particle, E/(NEp) 
and S/(Nks), as a function of the reduced temperature T '/Tp. 



2. Experimental measurement of equation of state 



Experimentally, there have been great efforts to measure the thermodynamics of strongly interacting Fermi gases of 
6 Li and 40 K atoms near a Feshbach resonance (26l - [32l |. Initial measurements have focused on trap averaged quantities 
[26M29l | . In the recent development, the bulk equ ation of state of a homogeneous Fermi gas becomes accessible (30l432l|. 
following a theoretical proposal by Ho and Zhou |l22| . Here we focus on the measurements performed by Nascimbene 
et al. at ENS (30l | and by Ku et al. at MIT (32| . These two precise measurements allow a quantitative comparison 
with the virial expansion predictions. 

In the ENS experiment, the local pressure P(fj,(z), T) of the trapped gas was directly probed using in-situ im ages of 
the doubly-integrated density profiles along the long z-axis (see the theoretical proposal by Ho and Zhou, ref. [122|). 
The temperature was determined by using a new thermometry approach employing a 7 Li impurity. The chemical 
potential could also be determined using the local density approximation, with fi(z) ~ [i — Vr (z) and the central 
chemical potential /1 being determined appropriately. By introducing a universal /i-function |l56l | 

h M = P ^ r ) (94) 

experimentalists were able to determine h(z) with very low noise. Here, P(/i,T) is the interacting pressure and 
P (1) (/^,T) is the pressure of an ideal two-componen t Fer mi gas. All the other thermodynamic quantities may then be 
derived from the universal h- function, i.e., see ref. |l23| . 
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In the MIT experiment (32|, instead of the pressure, the density equation of state p(p,T) — dP(p,T)/dp is 
measured. Owing to the perfect cylindrical symmetry of the trapping potential, the 3D density p{p, T) can be 
reconstructed from the measured column density, i.e., P2d(x,z) — j dyp(x,y, z), by using an inverse Abel transform 
(3^ . The local pressure and isothermal compressibility k = \p 2 dpldp\~ x can then be calculated from the density (32|. 
The crucial advantage of the MIT experiment is that the temperature T and the chemical potential p can be replaced 
by the pressure and compressibility. Thus, the notoriously difficult thermometry of a strongly interacting Fermi gas 
may not be required. 



3. Qualitative comparison between theory and experiment 




Figure 12: (color online) The second-order virial prediction for the ratio £i„t/ekin of a Fermi gas of 6 Li atoms near the 
Feshbach resonance. The dashed, solid, and dotted lines are for T/Tf — 1.2, 0.6, and 0.4, respectively. The symbols show 
the experimental data from the ENS group [26| . measured at T = 3.5 pK and T/Tf = 0.6. The squares and circles are 
respectively the results obtained by approaching the Feshbach resonace (Bo) from the BEC side and BCS side. From ref. [7(J 
with permission; copyright (2004) by APS. 



Before quantiative comparing the virial theory with the latest thermodynamics measurements, we mention briefly 
the first application of virial expansion in ultracold atomic Fermi gases, reported by Ho and Mueller in 2004 [7(3|. This 
elegant application gave a very good qualitative explanation for the measured interaction energy at ENS in 2003. 

Experimentally, in a Fermi gas of 6 Li atoms the Feshbach magnetic field was swept across the resonance from either 
the positive (BEC) or negative (BCS) scattering length side. The interaction energy of the near-resonance Fermi gas 
was then recorded at different fields. As shown in Fig. [T2] by symbols, crossing the resonance from the BCS side 
(the scenario A), the interaction energy remains negative and continuous across the resonance, while approaching 
the resonance from the opposite BEC side (the scenario B), £„( is positive but drops to a negative value near the 
resonance. 

By using the virial expansion theory to the second order, Ho and Mueller showed conclusively that the different 
interaction energy is a result of the different initial state. In the scenario A, the system is alway in the ground state 
with strong attractions, while in the scenario B, the system is initially in the metastable excited branch, where the 
interaction between two fermions is repulsive. For the detailed discussion, see Fig. O To be concrete, to the second 
order of virial expansion, the kinetic energy and interaction energy can be written as [70l |. 



3pk B T ( ^ p\\ B 

2 7/2 



(95) 



and 



tint 



3pk B T / 3 v 
— o — (AibJ 



1 dAb 2 
3 dT 



(96) 



The second-order virial coefficient Ab 2 can be calculated using the Beth-Uhlenbeck formalism Eq. (|3"Tj) and the usual 
s-wave phase shift. In the metastable excited branch, the contribution of the bound state to A62 should be removed. 
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The ratio tintl^kim predicted by Eqs. (|95[) and (|96[) . is compared with the experimental data in Fig. 1121 The virial 
prediction at T/Tp = 1.2 agree well with the experimental results, which were measured at T/Tp = 0.6. The difference 
in temperature is understandable, since the fugacity at T/Tp = 0.6 is already larger than 1 and the agreement must 
be affected by the higher order terms in the virial expansion. 



Quantitative comparison: Homogeneous system 
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Figure 13: (color online) Virial expansion prediction of the universal function h(z) up to the second (dashed line), third (solid 
line) and fourth order (thin solid line with error bar), compared with the experimental data (empty squares with error bars). 
Here, Ab 4 = —0.016 ± 0.004. Adapted from ref. I7p w ith inclusion of experimental results. The experimental data of the ENS 
group and of the MIT group are taken from refs. [30| and [3^ |. respectively. 



We now turn to the quantitative comparison. At high temperature, by using the virial thermodynamic potential 
Eq. (|92[) . the universal function h(z) can be written as, 



h{z) = 1 + 



Ab 2 z 2 + A6 3 2 2 



(2/y/F) / °° iVa In(l + «e-*)dt' 



(97) 



In Fig. 1131 we compare the virial expansion prediction and the experimental data for the universal function h(z). 
The virial results are calculated by using Eq. (|fJT|) . with inclusion virial coefficients Ab n up to Ab 2 = l/v2 (Virial2), 
Ab 3 ~ -0.35510298 (Virial3), and A6 4 = -0.016 ±0.004 (Virial4). At small fugacity {z < 0.7), the experimental data 
agrees excellently well with the virial prediction. Using A63 and A64 as independent fitting parameters, experimentally 
it was determined that Ab 3 , expt = -0.35 ± 0.02 and Ab^ expt = 0.096 ± 0.015 [30| by the ENS group. The latest 
measurement at MIT also showed an excellent agreement with virial expansion and reported A64 exp t = 0.096 ± 0.010 
[33 | . Thus, while the theoretical A63 has been confirmed unambiguously by the experiments, the theoretical prediction 
for the fourth virial coefficient A64 = —0.016 ±0.004 contradicts with the experimental observation. This discrepancy 
remains to be resolved. A possible reason is the uncertainty of the Feshbach resonance position, which is about 1.5 
G for 6 Li atoms [H]. As we discussed earlier, this uncertainty will leads to 1% relative error to the second virial 
coefficient. When this systematic error passes to the small fourth virial coefficient, the experimental determination of 
A64 may become unreliable. 

Let us turn to the other thermodynamic quantities such as energy and entropy. We report in Fig ll4l the temperature 
dependence of energy and entropy of a unitary Fermi gas in homogeneous space. The solid line and dashed line are 
the predictions of virial expansion up to the third-order and second-order, respectively. For comparison, we also show 
the ideal gas result by the thin dot-dashed line and the experimental results by symbols. We observe that the virial 
expansion is valid down to the degenerate temperature Tp , where the prediction up to the second-order or third-order 
expansion does not differ largely. The experimental data lie between the two virial expansion predictions, but clearly 
agree much better with the third-order expansion, as anticipated. 



27 




Figure 14: (color online) Energy per particle E/(NEf) and the entropy per particle S/(NEf) as a function of reduce tem- 
perature T/Tf for a homogeneous Fermi gas in the unitary limit. The predictions of virial expansion up to the second- and 
third-order are shown by dashed line and solid line, respectively. For comparison, we plot the ideal gas result by the dot-dashed 
line. We show also the experimental data measured at ENS [30] and MIT [32l |. which agree extremely well with the prediction 
from virial expansion. Adapted from ref. [73l |: copyright (2010) by APS. 



5. Quantitative comparison: Trapped system 




Figure 15: (color online) Energy per particle E/(NEf) and the entropy per particle S/(NEf) as a function of reduced 
temperature T/Tf for a trapped Fermi gas in the unitary limit. The predictions of quantum virial expansion up to the second- 
and third-order are shown by solid line and dashed line, respectively. For comparison, we plot the ideal gas result by the 
dot-dashed line. The experimental data measured at ENS (30l | and MIT [3^ | are shown by empty squares and solid circles, 
respectively. Adapted from ref. [73| : copyright (2010) by APS. 



Experimentally, the thermodynamics of a harmonically trapped Fermi gas in the unitary limit can be determined as 
well, from the measured universal /i-function. For the details, see ref. (66j. We present in Fig. [I5]the high-temperature 
expansion predictions for energy and entropy, and compare them with the experimental measurement. We find a much 
broader applicability of virial expansion: it is now quantitatively applicable down to 0.5Tp, as confirmed by the precise 
experimental data at ENS [13] and at MIT [32| . This is largely due to the much reduced higher order virial coefficient 
in a harmonic traps, i.e., Ab n ^T — n _3 / 2 A6„. At large n, the reduction factor of n~ 3 / 2 is fairly significant, implying 
a better convergence of virial expansion and hence a much wider applicability. 



6. Reliability of virial expansion 



To better understand the reliability of virial expansion, we show in Fig. [T^the fugacity as a function of temperature, 
for a homogeneous or trapped Fermi gas in the unitary limit. These two curves are determined from the experimental 
universal function h(z) measured at ENS and MIT. By setting z = 1 as the criterion for qualitative reliability, we find 
that the virial expansion should be applicable at T > 0.7Tp for a homogeneous unitary Fermi gas and at T > 0.5Tp 
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Figure 16: (color online) Fugacity z as a function of reduce temperature T/Tf, for a homogeneous unitary Fermi gas (solid 
circles) and for a trapped unitary Fermi gas (empty squares). The results, calculated from the experimental universal function 
h(z) (3(il. l32|. have a relative error at a few percents. 



for a trapped unitary Fermi gas. As the typical experimental temperature for a unitary Fermi gas is about 0.5Tp, 
we thus demonstrate clearly the virial expansion method is a very useful tool for understanding the properties of a 
normal, strongly interacting Fermi gas. 



G. Virial equation of state for a spin-population imbalanced Fermi gas 

We consider so far the balanced Fermi gas with equal mass a nd sp in-populations. The virial expansion is applicable 
as well to the imbalanced systems with either unequal mass |l24| or spin-populations |l25) . Here we focus on the 
latter case with unequal spin-populations and use virial expansion to obtain the high-temperature spin susceptibility 
of a unitary Fermi gas. 

In the presence of spin imbalance, it is necessary to introduce two fugacities z-f- = cxp(/i-j-/fcsT) and z^ = 
cxp(^j,/fcsT), and to distinguish different spin-configurations. Quite generally, we may write the thermodynamic 
potential as, 

00 n 

n = -k B TQi 4~ Wfc, (98) 

where b n ^ is the n-th (imbalanced) virial coefficient contributed by the configuration with n — k spin-up fermions and 
k spin-down fermions. It is easy to see that the imbalanced virial coefficients satisfy the relation b n ,k — b n . n -k and 
J2k=o bn,k — b n - 

The calculation of b n .k is straightforward, following the standard definition of thermodynamic potential. We rewrite 
the grand partition function Z =Trexp[— (H — [i^N^ — /iiN^/ksT] in the form, 

00 n 

z = EE4 l ~S fe ^< (") 

n=0 fc=0 

where Q n ,k is the partition function of a cluster that contains n — k spin-up fermions and k spin-down fermions. It is 
apparent that due to the symmetry in spin configurations we have Q n ^k = Qn.n-k- The imbalanced cluster partition 
functions satisfy as well a sum rule X)fe=o Qn.k = Qn- By expanding the thermodynamic potential £1 — — k B T hi Z 
into powers of the two fugacities, the imbalanced virial coefficients can then be expressed in terms of the cluster 
partition function Q n ,k- 
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1. Virial expansion of an imbalanced Fermi gas up to the third order 

To be concrete, let us consider the imbalanced virial expansion up to the third order. To this order, we may write 
the grand partition function as Z = 1 + X\ + x 2 + X3, where 

xi = ZfQi.o + ^Qi.i, (100) 

X-2 = Z^Q-2S) + z^ziQ 2 ,i + zfQ 2 ,2, (101) 

and 

x 3 = z^Q 3 , + z^J-^a.i + z t zfQ 3 , 2 + zfQ 3 , 3 . (102) 

By introducing a symmetric cluster partition function Q s n = Q n ,o = Qn,n and using the properties of Q n ,k, it is 
easy to show that Q{ = Qi/2, Q 2 j — Q 2 - 2Q 2 , and Q 3il = Q 3 , 2 = Q 3 /2 - Q|. Using ln(l + x\ + x 2 + x 3 ) ~ 
(xi + x 2 + x 3 ) — (x\ + 2x\X 2 )/2 + xf/3, after some algebra we obtain b n k (k < n/2), 

h,o = 1/2, (103) 

b 2 ,o = Q2/Q1 ~ Q1/8, (104) 

b 2 ,i = Q2/Q1 ~ 2Q S 2 /Q 1 - Qi/4, (105) 

b 3 ,o^Q 3 /Qi-Q S 2/2 + Qj/24, (106) 

and 

63,1 = Qs/(2Qi) - Q3/Q1 - Q2/2 + Q s 2 /2 + Ql/8. (107) 

The virial coefficients with k > n/2 can be obtained directly since b„_k = bn,n— k- 

As before, it is convenient to consider the interaction effect on the virial coefficients or the differences such as 
AQ n = Q n — Qrt , Ab n = b n — bn \ and Ab n ,k = b n ^ — b^L Here, the superscript "1" denotes an ideal, non- 
interacting system with the same fugacities and the operator "A" removes the non- interacting contribution. It is 
clear that the symmetric cluster partition function Q s n is not affected by interactions since the interatomic interaction 
occurs only between fermions with unlike spins. Thus, we have A62.0 = A63.0 = 0, A&2,i = A(&2 — 262,0) = A&2 and 
A&3.1 = A{b 3 /2 — 62,0) = A&3/2. Accordingly, we may rewrite the thermodynamic potential into the form (up to the 
third order), 



n = 0^ - k B TQ l 



z r z i Ab 2 + -A63 



(108) 



where fi^ 1 ) = f^ 1 ^/^) + fl^(pi) is the thermodynamic potential of a non-interacting Fermi gas. Hereafter, we 
consider the homogeneous case, in which Q\ = 2V/X 3 lB . 

With the virial expansion of thermodynamic potential Eq. (|108|) . we solve the standard thermodynamic relations 
iVf = —dQ/dfj^ and iVj, = —dQ/dp,^. for the two fugacities z\ and z±, at a given reduced temperature r = T/Tp and 
a given spin imbalance P = (N^ — N±)/N. Here, Tp = fi, 2 (37r 2 p) 2 / 3 /(2m)//cB is the Fermi temperature. It is easy to 
show that we can define a dimensionless number density p = pX 3 /2 = 4/ (By^T 3 / 2 ), p^ = (1 + P)p, and p_i = (1 — P)p. 
We then rewrite the number equations into dimensionless forms, 

p t = (z t ) + z t z i 2A6 2 + (2.^ + z t zf) A63, (109) 
Pi = P il Hzi) + z t z l 2Ab 2 + {z^z l + 2z t zl)Ab 3l (110) 

where ,5« (z) = (2/v^F) Vt[ze~ 1 / (1 + ze-*)]dt. We can obtain the two fugacities by solving the coupled number 
equations. 
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2. Virial expansion of spin susceptibility and compressibility 

We now calculate the spin susceptibility xs = (dSp/dSp) in the balanced limit of P = 0. For this purpose, we 
determine the two by two susceptibility matrix S = {dp a / dp a i) to the third order of fugacity. For a homogeneous 
unitary Fermi gas, using the number equation we find that, 



S(P = 0) 



1 



k B TX 3 



A B 
B A 



[III) 



where 



\/TT 



(1 + ze~ t ) 



rdt + 2z 2 Ab 2 + 5z 3 Ab 3 



(112) 



and 



B = 2z 2 Ab 2 +4z 3 A6 3 . 



(113) 



The spin susceptibility xs = 2(^4 — B) / (kBTX 3 dB ) and compressibility k = 2(A + B) / (kBTX 3 dB ) are then given by, 



k B TX 3 



s/tze-* 
7T J fl + ze"*) 





-dt + z 3 Ab 3 



(114) 



and 



k B TX 3 



dB 



y/ize-* 



J (l + ze- 1 ) 
o 



Trdt + 4z 2 Ab? + 9z 3 Ab 3 



(115) 



respectively. In the expressions, we see clearly the effect of interactions. For the spin susceptibility, it appears in the 
third order of fugacity only. A high-temperature measurement of spin susceptibility therefore could be a sensitive 
way to measure accurately the third virial coefficient. We note that, for an ideal Fermi gas the spin susceptibility and 
compressibility are equal. 
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Figure 17: (color online). Spin susceptibility and compressibility of a homogeneous Fermi gas in the unitary limit, normalized 
by the T = ideal Fermi gas susceptibility value \o — K o = 3p/(2Ef)- Here p is t he density and Ef the Fermi energy. The 
experimental compressibility data are taken from ref. [32l |. Adapted from refs. [l25l | and [32l |. 



Fig. [T7]reports the numerical result of Eqs. (| 1 14[) and (|115p for a homogeneous Fermi gas in the unitary limit, where 
the fugacity z is solved consistently to the third order expansion in the number equation. The spin susceptibility 
and compressibility are smaller and larger than that of an ideal, non-interacting Fermi gas, respectively, as expected. 
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Experimentally, the spin susceptibility at finite temperatures is related to the measurement of the thermal spin 
fluctuations: 

A (N-f - Ni) 2 = xs 

N p K ' 

A shot noise measure ment of the spin fluctuations therefore could be used as a sensitive thermometry for strongly 
interacting Fermi gases |l26j . provided that the spin susceptibility is known. Now, this seems to be accessible, since the 
shot noise me asuremen ts of the density fluctuations in a weakly interacting Fermi gas have already been demonstrated 
very recently [L^[i"28l|. We find that the spin susceptibility is strongly suppressed by interactions with respect to the 
ideal Fermi gas result, even well above the degenerate temperature Tp. At T = Tp, the reduction is about 40%. For 
the compressibility, the virial prediction agrees very well with the latest experimental measurement at T > Tp [32], 
as anticipated. 



III. VIRIAL EXPANSION OF TAN'S CONTACT 



In this section, we show that the important many-body parameter - Tan's contact - can be virial expanded in terms 
of the so-called contact coefficients [77j]. By using few-particle solutions, we determine the second and third contact 
coefficients. For a trapped Fermi gas in the unitary limit, we find that the virial prediction at T > 0.5Tf agrees very 
well with the recent experimental measurements performed at Swinburne University of T echnology |l07l Il08| . The 
first virial expansion calculation of Tan's contact was given by Yu, Bruun, and Baym |l29( . 



A. Virial expansion of Tan's contact 



The virial expansion of the contact follows directly from an alternative representation of Tan's adiabatic sweep 
theorem in the grand- canonical ensemble, 



on 



di-aj 1 ) 



Airm 



-I. 



(117) 



(118) 



This is simply because the adiabatic sweep theorem implies the first law of thermodynamics, 

AE = h 2 l/{4.Tzm)A{-a~ l ) + TAS + (lAN , 
which can alternatively be written as 

An = ft 2 I/(47rm)A(-a ; : 1 ) - SAT - NAfi . (119) 
Therefore, using virial expansion for £1 we immediately obtain a quantum virial expansion for the contact: 

4nmkBT\dB 



1 = 



h c 2 z* + 



C n z' 1 + 



(120) 



where we have defined the dimensionless contact coefficient, c n = dAb n /d(XdB/o-s)- For a homogeneous system, we 
shall use the contact intensity, C ~ I/V . 

In general, the contact coefficient should be a function of Ads/a s and hence is temperature dependent. In the 
unitarity limit where A<2s/a s =0, however, we anticipate a constant, universal contact coefficient, similar to the 
universal virial coefficient Ab n WL 72| . This is a manifestation of fermionic universality, shared by all systems of 



strongly interacting fermions [lOL lll| . 



B. Universal relation between homogeneous and trapped contact coefficients 

In exact analogy with the virial coefficient, fermionic universality leads to a very simple relation between the trapped 
and homogeneous contact coefficients at unitarity. Let us consider the contact of a harmonically trapped Fermi gas 
with the trapping potential Vr(x) = muJ T (x 2 + y 2 + z 2 )/2. In the thermodynamic limit of uip — > 0, as before we may 
use the local density approximation and neglect the discrete energy levels. The whole Fermi system is treated as many 
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cells with a local chemical potential /i(x) = fi — Vr(x) and a local fugacity z(x) = e M( x )/ fe B T = z exp[— Vr(x)/(fcsT)]. 
Due to the constant contact coefficients, the spatial integration in the total contact It = J dx[C(x)] can be easily 
performed. We find that, 

It = AlTmk f XdB Q^T [c,tz 2 + c 3 ,tz* + ...}, (121) 
where the trapped contact coefficient is given by a universal relation, 

Cn,T = (122) 

and Qi,t = 2(fcsT) 3 /(/iWT) 3 is the single-particle partition function in harmonic traps and in the local density 
approximation. 

In the following, using the known solution of two- and three-fermion problems, we calculate the universal second 
and third contact coefficients, in both homogeneous and trapped configurations. 

C. Second contact coefficient 

The second contact coefficient of a homogeneous interacting Fermi gas can be obtained from the well-known Beth- 
Uhlenbeck formalism for the second virial coefficient. In the vicinity of the unitary limit, we have Ab2(a s < 0) ~ 
l/y/2 + Ads/(7ra s ), giving rise to a homogeneous contact coefficient, 

c 2 = -. (123) 

To calculate the trapped second contact coefficient, we consider the second virial sufficient in an isotropic harmonic 
trap, which is given by Eq. (J52J), Ab 2 , T = (1/2) Y, Vn [e^^ +3 / 2 ^ T - e -(2^ 1) +3/2)a r ] ) where q t = foj T /(/- B T) < 1 i s 
the reduced trapping frequency, v n satisfies the secular equation 2r(— v n )/T{— v n — 1/2) = d/a s , and d = y^KJJrrujJr) 
is the characteristic length scale of the harmonic trap. In the non-interacting limit, Vn = n (n = 0, 1, 2, ...), and in 
the unitary limit, v n = n — 1/2. It is easy to show that, 



d{\dB/a s ) 
Thus, we find that in the unitary limit, 



Ad B /a„=0 



27rAds 



*jz±±m. (i 2 4) 



r(n+ l/2) p _ (2 „ + i/2)a 

=o 

The sum over n can be exactly performed, leading to 



C2.T = ~z r lut E 1 , e ' (2 " +1/2) ^- ( 125 ) 

n—0 



l 

C2,T 



2ujt 



,+ujt 



1/2 

(126) 



2V2tt 



i-f + o(«4) 



The leading term in the above equation is universal, satisfying the universal relation Eq. (|122p . The second term 
(oc d> T ) is non-universal and is caused by the length scale of the harmonic trap It represents the finite-size 

correction to the local density approximation that we have adopted above. 

D. Third contact coefficient 

The determination of the third contact coefficient is more cumbersome. As in the calculation of the third virial 
coefficient, we can determine firstly the trapped contact coefficient 03^, and then to use the universal relations at low 
trap frequency to obtain the homogeneous result, C3 = 3-\/3cs t T- 

An estimate of C3,t can already be obtained by the known results of Ab^^T as a function of the coupling constant 
l/kpa-s at different temperatures T/Tp and lot ~ 0.15, as shown in Fig. 8. This is simply because, 



1 dAb 3 . T / T dAb^ T 
° 3 T = k F X dB d{l/k F a s ) ~ V 4*T F d(l/k F a s y (12?) 
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We find that the coefficient C3.T at resonance is indeed nearly temperature independent and estimate from the 
slope of Ab^ t T that, c^ t T (estimate) ~ —0.0265 at lot ~ 0.15. An accurate determination of c 3 x requires a 
systematic extrapolation to the limit of lot = 0. For this purpose, we calculate numerically the derivative 
C3,t(wt) = [fA^r 1 'd(XdB I ' a s)\x dB / O ,=o as a function of lot- Using the small ujt data, a numerical extrapola- 
tion to u>t = gives rise to the trapped third virial contact coefficient, c^t — —0.0271 ± 0.0002. Thus, we obtain 
immediately from the universal relation, Eq. (|122j) . the homogeneous third contact coefficient, 

c 3 = -0.1408 ±0.0010. (128) 

0.0 
-0.1 
-0.2 

o 

-0.3 
-0.4 



-0.5 

-2-10 1 2 

dB s 

Figure 18: (color online). Third contact coefficient as a func tion of the dimensionless parameter \<ib/cl b . The result is calculated 
using the third virial coefficient reported by Leyronas [lld |. The inset shows the second contact coefficient. 

Alternatively, we can determine the third contact coefficient by tak ing a numerical derivative of the third virial 
coefficient Ab^(XdB/o, s ), which was calculated recently by Leyronas [llOj, by using diagrammatic field theoretic 
method. As shown in Fig. [TSl in the unitary limit we find that C3 = —0.1399 ± 0.0001, in excellent agreement with 
the result Eq. (|128| from the exact three-particle solutions. 




E. Large-T contact: the homogeneous case 



We are now ready to calculate the universal contact in the high temperature regime. For a homogeneous Fermi 
system, the single-particle partition function Q\ = 2V/X 3 iB and the dimensionless contact I/(Nkp) is given by, 

aIf^ 3 * 2 ^) 2 ^ 2 ^- 1 ' <129) 

Here N = pV is the total number of atoms with the homogeneous density p. 
The fugacity z is determined by the number equation [73|, 

p = (z) + [2A6 2 z 2 + 3A6 3 z 3 + •••], (130) 

where we have defined a dimensionless density p = p\ dB /2 = [4/ (3\/tt)}(Tf /T) 3 ' 2 and the density of a non-interacting 
Fermi gas as p^{z) = (2/0F) / °° dtVt/(l + z~V) . 

In practice, for a given fugacity, we calculate the dimensionless density using Eq. (|130[) and hence the reduced 
temperature T/Tp. The dimensionless contact is then obtained from Eq. (|129|) . as a function of T/Tp or the inverse 
fugacity z . 

Fig. [19] reports the temperature (main figure) or fugacity (inset) dependence of the homogeneous contact in the 
unitarity limit, calculated by virial expanding to the second order (dashed line) or third order (solid line). The close 
agreement between the second and third predictions strongly indicates that the virial expansion works quantitatively 
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Figure 19: (color online) Universal contact of a homogeneous Fermi gas in the unitary limit at high temperatures, as predicted 
by the virial expansion method up to the second order (dashed line) and the third order (solid line). Dashed vertical line 
indicates the Fermi degenerate temperature Tf ■ The inset shows the contact as a function of the inverse fugacity. From ref. 




well down to the Fermi degenerate temperature Tp, as indicated by the vertical dashed line. At sufficient high 
temperatures, where 

z^=[4/(3v^)](T f /T) 3 / 2 , (131) 
the leading temperature dependence of the contact is given by, 

SHS^ 1 " -"(£)"■ (132) 

as predicted by Yu and co-workers |l29| . We note however that the pre- factor there is smaller by a factor of An 2 , due 
to a different definition for the contact. 



F. Large-T contact: the trapped case 



For a trapped Fermi gas at unitarity, the dimensionless contact can be written as, 




(133) 



where c 2 , T = l/(2y/2n) and c 3 , T = -0.02692 ± 0.00002. The number equation takes the form 

p T = p { t ] (z) + [2Ab 2 . T z 2 + 3A6 3 , T z 3 + ■■■], (134) 

where pr = (N/2)(Hujt) 3 /(ksT) 3 = (Tp /T) 3 /6 and the density of a non-interacting trapped Fermi gas Pp\z) = 
(1/2) J °° dtt 2 /(l + z~ 1 e t ). hi analogy with the homogeneous case, for a given fugacity we determine the reduced 
temperature T/Tp from the number equation (|134|) and then calculate the trapped contact using Eq. (|1 33[) . 

Fig. [20]presents the virial expansion prediction for the trapped universal contact, expanding up to the second order 
(dashed line) or third order (solid line). Amazingly, because of the factor of ?i~ 3//2 reduction for the n-th contact 
coefficient in harmonic traps, the convergence of the expansion is much improved. The expansion now seems to be 
quantitatively reliable down to 0.5Tp. The asymptotic behavior of the contact at very high temperatures can be 
determined by setting 



z — Pt = {Tp/Tf/Q . 



(135) 
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Figure 20: (color online) Universal contact of a trapped Fermi gas in the unitary limit at high temperatures, obtained by 
expanding the virial series to the second order (das hed line) and the third order (solid line). The inset shows the contact as a 
function of the inverse fugacity. Adapted from ref. [107| | ; copyright (2011) by APS. 



We find that. 



Nk, 



(T > T F ) 



(i 



T 
TV 



-5/2 



(136) 



Thus, the contact in harmonic traps decays at high temperatures much faster than in homogeneous space, due to the 
reduction of the peak density at the trap center at high temperatures. 

The finite-temperature contact of a trapped Fermi gas in the unitary limit was recently measured at Swinburne 
University of Tec hnol ogy. Using the structure factor Tan relation Eq. (fTf5| , the contact was extracted from the static 
structure factor |l07| . which has been measured by two-photon Bragg spectroscopy. In Fig. I20[ the experimental 
result was shown in solid circles. At T > 0.5Tp, the data agree well with the virial prediction. 



IV. VIRIAL EXPANSION OF DYNAMIC STRUCTURE FACTOR 

So far we consider the virial expansion of static properties of a strongly correlated Fermi system. In the following, 
we show that dynamic properties can be studied as well using virial expansion. This issue is less explored in the 
literature. In this section, we consider the dynamic density response of a strongly correlated Fermi system [75| . 

A. Dynamic structure factor 



The dynamic density response is characterized by the so-called dynamic structure fact or ( DSF), which gives the 
linear response of the many-body system to an excitation process that couples to density |95| . For ultracold atomic 
gases, it can be conveniently measured by two-photon Bragg spectroscopy using two laser beams (48| . Theoretically, it 
is difficult to predict DSF in the strongly interacting regime. Traditional tools such as the perturbative random-phase 
approximation (RPA) theory are in principle reliable in the weakly interacting limit |l3lMl3 4| . 

The DS F S( q, io) is the Fourier transform of the density-density correlation functions at two different space-time 
points [95l Il3dl |. For a balanced atomic Fermi gas with equal spin populations N/2 (referred to as spin-up, a =f, and 
spin-down, a =1), Syf(q,u) = Sn(q, lu) and S|j,(q, w) = 5j.^(q, w), each of which is defined by, 



S aa ,(q,Lu)=Q- 1 Y,e~ En ' /kBT '{n\5 Pr7 (q)\n')U' Spl(q) n) S {tvjj - E nn ,) 



(137) 



Here \n) and E nn i = E n — E n i are, respectively, the eigenstate and eigenvalue of the many-body system, while 
Q = J2 n exp{—E n /kBT) is the partition function. The density operator 5/3 CT (q) = ^2 ia e~ ?q Xi is the Fourier transform 
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of the atomic density operator Sp a (r) for spin- g- ato ms. The total DSF is given by S(q, u>) = 2[S^(q, uj) + S^i(q, u>)]. 
The DSF satisfies two remarkable /-sum rules |135| . 



(138) 



+ °° hq 2 

S(q, u>)ujduj = N 

.00 2m 



and 



S^(q,u)ujduj = 0, 



(139) 



which hold irrespective of interactions and temperatures. 

According to the finite-temperature quantum field theory |l30l |. it is convenient to calculate DSF from dynamic 
susceptibility, Xaa' (Qj t ) = — {T T f>a (q^ r ) Pa' (q 5 0)), where r is an imaginary time in the interval < r < j3 = 1/kgT. 
The Fourier component Xaa 1 (q, *w n ) at discrete Matsubara imaginary frequencies iuj n = i2nirkBT (n = 0, ±1, ...) gives 
directly the DSF, after taking analytic continuation and using the fluctuation-dissipation theorem: 



S aa > (q,w) = - 



ImXgg' (q; -> a; + ^0+) 

TfM — e -Hu/k B T\ 



(140) 



The frequency integral of the DSF defines the so-called static structure factor (SSF). For different spin components, 
we have, 



2 r+00 

Saa' (q) = -T7 / Sao-' (q,w) dw. 

^ —00 



(141) 



The total SSF is given by, S (q) = (1/iV) JJ^ dwS' (q,w) = 5ft (q) + %(q). As we mentioned earlier, the SSF is 
related to the two-body pair correlation function g aa / (r) (95| through a Fourier transform. 

Experimentally, the DSF is measured by inelastic scattering experiments of two- photon Bragg spectroscopy [48[. 
The atoms are exposed to two laser beams with differences in wave-vector and frequency. In a two-photon scattering 
event, atoms absorb a photon from one of the beams and emit a photo into the other. Therefore, the difference in 
the wave-vectors of the beams defines the momentum transfer fiq, while the frequency difference defines the energy 
transfer hui. In the regime of large transferred momentum, which is exactly the case in current experiments for the 
crossover Fermi gas [48| . the single-particle response is dominant and peaks at the quasi-elastic resonance frequency 
uj res = hq 2 / (2M), where M is the mass of the elementary constituents of the system. Therefore, we may anticipate 
that the Bragg response peaks at ojr = hq 2 /(2m) in the BCS limit and peaks at ujR, m oi = hq 2 /(Am) = ujr/2 in the 
BEC limit, since the underlying particles are respectively free atoms (M = m) and molecules (M = 2m). 
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Figure 21: (color online) (a) Measured dynamic structure factor of a harmonically trapped Fermi gas in the BEC-BCS crossover 
at the lowest attainable temperature (T < O.ITf) and at a large transferred wave- vector q — The inset shows the static 

structure factor as a function of the dimensionless interaction parameter. Adopted from ref. with permission; copyright 
(2008) by APS. (b) Temperature dependence of dynamic structure factor of a trapped Fermi gas in the u nita ry limit, measured 
at q = 2.7&F. The inset shows the temperature dependence of static structure factor. Adapted from ref. |l07l |: copyright (2011) 
by APS. 
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In Fig. 1211 we summ arize the main experimental results for a harmonically trapped Fermi gas in the BEC-BCS 
crossover [48l, Il06l -[l08]. Fig. l2~Tk shows the DSF (main panel) and SSF (inset) at several dimensionless interaction 
strengths and at the lowest experimentally attainable temperature (i.e., T < O.lTp, where Tp is the Fermi temper- 
ature ) I48L wh ile Fig. 121b presents the temperature dependence of structure factors in the most interesting unitary 
limit |l07lll08| . As anticipated, in Fig. the DSF peaks at ujr/2 and uj r on the BEC side (i.e., l/(k F a s ) = +0.5) 
and on the BCS side {l/{kpa s ) = —0.8), respectively. In the unitary limit, where the statistics of the elementary 
excitations is not well defined, we observe a two-peak structure with responses from both molecules and free-atoms. 
As the temperature increases (Fig. 121b ). however, these two peaks merge. The resultant broad peak shifts eventually 
to ujr at high temperatures. 



B. Virial expansion of dynamic structure factor 

We construct first the virial expansion for the dynamic susceptibility Xaa' (x, x';r > 0), which is given by, 

Xaa ' = Tre -CH-M0/fc B T • ( 142 ) 

At high temperatures, Taylor-expanding in terms of the powers of small fugacity z = exp(/x//c_BT) <C 1 leads to 
Xaa' (x, x'; r) = {zX\ +z 2 X2 + - ■ ■ )/ (l + zQi + z 2 Q2 + - ■ ■ ) = zXi+z 2 (X2 — XiQi) + - ■ ■ , where we have introduced the 
cluster functions X n = — Tr n [e~ w / fcBT e' rW n cr (x)e _rW n (T ' (x')] and Q n =Tr„[e -W / fcBT ], with n denoting the number of 
particles in the cluster and Tr„ denoting the trace over n-particle states of proper symmetry. We shall refer to the above 
expansion as the virial expansion of dynamic susceptibilities, Xaa' (x, x'; r) = zxaa\i (x, r'; r)+z 2 Xaa' .2 (x, x'; r) + - • • , 
where, 

Xaa's (x,x';r) = Ai, 

X CTCT ', 2 (x,x';r) = X 2 -X 1 Q 1 , etc. (143) 
Accordingly, we shall write for the dynamic structure factors, 

Saa> (q,w) = zSaa',1 (q, oj) + Z 2 S a a',2 (q, H . (144) 



C. Trapped virial dynamic structure factor up to the second order 

The calculation of the n-th ex pan sion coefficient requires the knowledge of all solutions up to n-body, including both 
the eigenvalues and eigenstates |72l. [73l| . Here we aim to calculate the leading effect of interactions, which contribute 
to the 2nd-order expansion function [75) . For this purpose, it is convenient to define Axaa',2 = {Xaa 1 ,2} — {-^2}^ 
and /S.Saa> ,2 = {S a a' .2}^ ■ The notation {}^ means the contribution due to interactions inside the bracketed term, 
so that {A^}^ 7 ' = X2 — X^i where the superscript "1" in X^ denotes quantities for a noninteracting system. We 
note that the inclusion of the 3rd-order expansion function is straightforward, though involving more numerical effort. 

It is easy to see that Axaa',i = 0, according to the definition of notation {}^ . To calculate the 2nd-order 
expansion function for the dynamic susceptibility, Axaa',2 = — {Tr^i \e~' H l kBT e^h^ (x) e~ rn n a i (x')] }^\ we insert 
the identity J^q \Q) (Q\ = 1 an d take the trace over the state P, i.e., 

AX^',2 = [e- E rl k * T +^r-E Q ) ( p |^| Q) {Q | p) j W _ (M5) 

P.Q 

Here, P and Q are the two-atom eigenstates with energies Ep and Eq, respectively. Expressing the density operator 
in the first quantization: (x) = ^\ 5 (x — x^) and n± (x) = ^ . S(x — x^), it is straightforward to show that, 

Axaa',2 = £ {e- E ^ T +^- E ^ (x, x')} (/) , (146) 

where 

C^ Q = [ dx 2 dx' 2 [$£$ Q ] (x, x 2 ) [$^$ p ] (x', x' 2 ) (147) 
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and 



C^ Q = J d Xl dx 2 [$^ Q ] (x,x 2 ) [$£* P ] (x ljX ') 



(148) 



The dynamic structure factor can be obtained by taking the analytic continuation. This result is AS a a',2 ( x > x '; w ) = 
J2pq |^ ( w + Bp — Eq) e~P Ep C^®, (x,x')| . Applying a further Fourier transform with respect to r = x — x' and 
integrating over X = (x + x')/2, we obtain the response AS a(J >,2 (q 5 ui), 



AS^,2 =J2{H" + Ep- E Q ) e~ f3E -C P J (q)} 



(/) 



(149) 



P.Q 



where C%! (q) = / dxdxfe-^^-^C^ (x,x'). 

The calculation of C P ® (x, x') or (q) is straightforward but tedious, by using the two-atom solutions in an 
isotropic harmonic trap mui 2 ^x 2 /2. We refer to ref. [75fl for details. The final result is given by, 



V p2q2 

where B = (k B T) 5 / 2 / (qh 4 ^), and 

A^ 2 = (-1)'( 1 -^')(2Z + 1) 



1 



p2(}2 



(150) 



drr2jl (f~) ^ ^ n «'« ^ 



(151) 



In Eq. (|15ip . we specify p2 = {n p l p } and q2 = {n q l q }, and I = max{Z p ,^ g } for the two-atom relative radial wave 
functions 4>{r) with energy e. We require that either l p or l q should be zero (i.e., minjZp,^} = 0), otherwise A p 2 q2 
will be cancelled exactly by the non-interacting terms. 



Together with the non-interacting DSF S^, , we calculate directly the interacting structure factor, 

S aa . (q, u) = Sfj, (q, u) + Z 2 AS <7<7 ' 1 2, 



(152) 



once the fugacity z is determined by the virial expansion for equation of state. 



1. Comparison of theory with the Swinburne experiment 

Considerable insight into the dynamic structure factor of a strongly correlated Fermi gas can already be seen 
from Eq. (|150p . in which the spectrum is peaked roughly at Wr iOTO / = lor/2, the resonant frequency for molecules. 
Therefore, the peak is related to the response of molecules with mass M = 2m. Eq. (|150p shows clearly how the 
molecular response develops with the modified two-fermion energies and wave functions as the interaction strength 
increases. In the BCS limit where AS arT i 2 is small, the response is determined by the non- interacting background 
that peaks at cup. In the extreme BEC limit (a — > + ), however, AS a(J i t 2 dominates. The sum in AS aa ' t 2 is 
exhausted by the (lowest) tightly bound state 4> re i(r) ~ \[2j 1 a s e~ r l aa with energy e re i ~ Ep = —h 2 /(ma 2 ). The 
chemical potential of molecules is given by fi m = 2/i — Ep. Therefore, the DSF of fermions takes the form, 



- z m B\/—exp 



J (7(7 



(W - UJR,mol) 
4UR,molKB-l 



where z m — e tJ ' m ^ kBT is the molecular fugacity. This peaks at the molecular resonant energy. As anticipated, 
Eq. (|153p is exactly the leading virial expansion term in the DSF of non-interacting molecules. It is clear that 
Sft(q, lj) ~ Sf4.(q, w) in the BEC limit, since the spin structure in a single molecule can no longer be resolved. 

To understand the intermediate regime, in Fig. !2iZb we report numerical results for the DSF as the interaction 
strength increases from the BCS to BEC regimes at T = 0.5Tp [75j]. The temperature dependence of the DSF in 
the unitary limit is shown in Fig. [22b [75]. In a trapped gas with total number of fermions N, we use the zero 
temperature Thomas-Fermi wave vector k p = (24iV) 1 / 6 /a/ IO and temperature Tp = (SNy^hux/kp as characteristic 
units. In accord with the experiment (48l . Il07j . we take a large transferred momentum of q = 3kp. At T = 0.5Tp, A 
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Figure 22: (color online) (a) Evolution of dynamic structure factor of a trapped Fermi gas in the BEC-BCS crossover with 
increasing interaction strength l/(kFa a ) at T = 0.5Tf. (b) Temperature dependence of dynamic structure factor of a trapped 
unitary Fermi gas. The dark circles indicate the peak position of spectra. The transferred wave- vector is q = 3/cf- Adapted 
from ref. [H; copyright (2010) by APS. 
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Figure 23: (color online) Comparison between theory and experiment for the dynamic structure fa ctor of a trapped unitary 
Fermi gas at finite temperatures. Here, the transferred wave-vector is q = 2.7kp. Adapted from ref. [l08l ] : copyright (2011) by 
NJP. 



smooth transition from atomic to molecular responses is evident as the interaction parameter l/(kpa s ) increases, in 
qualitative agreement with the experimental observation (c.f. Fig. l2"Tk). In the unitary limit, the peak of total DSF 
shifts towards the molecular recoil frequency, as indicated by the dark circles. This red-shift is again in qualitative 
agreement with experiment (c.f. Fig. 121b ). 

For a close comparison, we plot in Fig. [23] the virial expansion prediction and experimental data for the DSF at 
several temperatures in the unitary limit. The theory is in very go od agreem ent with experimental data at high 
temperatures (see, for example, the case of T = 0.69Tf in Fig. I2"3"k ') |l07l Il08| . where the fugacity z is less than 1. 
Towards low temperatures, the agreement becomes worse. However, the virial expansion does capture the qualitative 
feature of the DSF, for temperature down to the onset of superfluid transition, T c ~ 0.2TV- 



D. Homogeneous virial dynamic structure factor up to the second order 

Let us now consider the expansion functions of a homogeneous Fermi gas in the unitarity limit. This can be extracted 
from the trapped expansion function because of fermionic universality in the unitary limit. As the scattering length 
diverges, all microscopic scales of the interaction are absent (Io| . For this few-body problem, the only energy scale is 
fcgT and length scale is the thermal de Broglie wavelength XdB- Dimensional analysis leads to, 



v 

A5 (T(T ',„(q,a;,T) = 3 AS aa >, n (q,&), 



(154) 
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where V is the volume, q = [S 2 q 2 /(2mfcsl 1 )] 1 / 2 , Co = fuo / '(fc^T), and A5 ffo ' in is a dimensionless expansion function. 
The temperature T is now implicit in the variables q and Co. This universal form implies a simple relation between the 
trapped and homogeneous expansion function. In a shallow harmonic trap, Vr(x.) = mio^ix 2 +y 2 + z 2 )/2, where lot 
0, the system may be viewed as a collection of many cells with a local chemical potential /i(r) = \i — Vt{y) and fugacity 
z(r) = zexp[-V T (r)/k B T], so that the trapped DSF is given by AS aa >,T (q,w,T) = / dr[AS cc > (q,w,T,r) /V]. Owing 
to the universal q- and w-dependence in the expansion functions, the spatial integration can be easily performed, giving 
rise to 



AS aa ., n {q,u) = r?' 2 



(k B Tf 



(q.w.T). 



(155) 



The (non-universal) correction to the above local density approximation is at the order of 0[(hcoT) 2 / (ksT) 2 ]. Eq. 
(|155D is vitally important because the calculation of expansion functions in harmonic traps is much easier than in free 
space. 
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Figure 24: (color online) Universal second order expansion function of DSF at q = 1/3, 1, and 3. The inset shows the rapid 
convergence of A5 , 2(q, to) at small hwr/kBT (thick lines) and, A5-|-f,2 (thin solid line) and A5^2 (thin dashed line) at q = 1. 
From ref. H3|. 




Fig. [M] reports the homogeneous expansion function AS2 = 2[A5^,2 + AS-j-j.^] at three different momenta, using 
AScrcr',2,T in Ref. [75[ as the input. One observes a quasielastic peak at Co = q 2 /2 or to = ftq 2 /(4m), as a result of 
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the formation of fermionic pairs. In Fig. 1251 we show the total homogeneous dynamic structure factor at q = 5/cf, 
calculated up to the second order. 



1. The f-sum rules 



We may derive sum rules that constrain the expansion functions, using the well-known /-sum rules satisfied by 
DSF. Using the virial expansion of the total number of fermions N, we shall have /-sum relations 

u}AS^. n (q, ui)dui — nq 2 Ab n (156) 

J — oo 

and 

" + 00 

u)AS^, n (q, &)du) = 0, (157) 

o 

which hold for arbitrary transferred momentum. 



2. Virial and contact coefficients from the large-q expansion functions 

At large momentum, the spin-antiparallel static structure factor satisfies the structure factor Tan relation Eq. ()16l) , 
/ <S^(q, U), T)du) ~ 1/ (8hq). By virial expanding both sides of the equation, we find that, 

r+°° ~ n 3/2 
AS n , n (q » 1) = / AS n , n (q,u>)du> = (158) 

On the other hand, in the same lim it of large momentum, the spin-parallel static structure factor is nearly unity so 
that / S n (q,u,T)du ~ AT/ (2ft) 0,[lO6|. This leads to 

/+oo 
AS ntn (q,u))dQ = nAb n . (159) 
-oo 

For the second expansion function, AS a(7 ^2, we have checked numerically that all the above mentioned sum rules are 
strictly satisfied. 



V. VIRIAL EXPANSION OF SINGLE-PARTICLE SPECTRAL FUNCTION 



In this section, we present the virial expansion of single-particle spectral function, a quantity that plays a key role 
in understanding the nature of pairing in strongly correlated Fermi gases. It has been argued that there might be a 
small window for pseudogap - the precursor of fermionic pairing in the normal state above the superfluid transition 
temperature - in analogy with high-T c superconductors [52[. However , its unambiguous identification is still unde r 
debate. Some of strong-coupling theories predict a pseudogap |l36l - ll40l |. while some other s claim no such effects |l4l| . 
Ab-initio quantum Monte Carlo simulations of the spectral function have been performed |l42Lll43ij . but the accu racy 
is yet to be improved. To date, the experimental measurements, throu gh the momentum- resolved rf spectroscopy |46| . 
were not conclusive, though there is a weak indication of pseudogap [47J. Here, we show that one can use the virial 
expansion up to the second order to qualitatively understand the experimental results (76| . Further improvements of 
virial expansion might be useful to solve the delicate pseudogap problem. 



A. Virial expansion of single-particle spectral function 



To virial expand the single- par ticle spectral function, let us consider the related finite-temperature Green function 
at different space-time points [76|], 



Geo' ( x > x '; r ) 



Tr 



^ Te -CH-fj,X)/k B T 



(160) 
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where at finite temperatures we are working with an imaginary time r in the interval < r < f3 = l/fc^T. At high 
temperatures, both numerator and denominator may be expanded into the powers of z <C 1, leading to G CT(T < (r, r'; r) = 
(X +zX 1 + --- )/(l + z Q 1 + ...)=X + z(X 1 - XqQx) + ■■■, where X n = - Tr„[e-«/ fc « T * (T (x,r) *+,(x')] is the 
expansion function and Q n =Tr n [e~ w / fcsT ] is the cluster partition function. The above expansion is to be referred to 
as the virial expansion of Green function, G aa ' (x, x'; r) = G<t<t',o (x, x'; t) + zG aa > ,i (x, x'; r) H , where, 

(x,x';r) = ^o, (w,i (x,x';r) = Xi - X Qi, etc. (161) 

We then take the Fourier transformation with respect to t, to obtain G aa i (x, x'; itu n ). The experimentally measured 
spectral function A (k, w) can be calculated from the finite-temperature Green function via analytic continuation, 

A aa i (x,x'; w) = --ImG ffff - (x,x';i^« -> w + i0 + ) . (162) 

7T 

A final Fourier transform on x — x' leads to A aa -i (k, uj) , as measured experimentally. For a normal, balanced Fermi 
gas, Af-f = = A(k, w) and _A-|-j, = 0. In accord with the virial expansion of Green function, we may write the 
spectral function, 

A (k, u) = A (k, u) + zA x (k, u) + ■ ■ ■ . (163) 

The calculation of the n,-th expansion function G n (x, x'; r) or A n (k, u>) requires the knowledge of solutions up to the 
n-body problem, including both energy levels and wavefunctions. 

As before, in the calculations of the Green function or spectral function, it is convenient to separate out the 
contribution arising from interactions. To this aim, for any physical quantity Q we may write Q = {Q}' 7 ' + Q^, 
where the superscript "1" in denotes the part of a non-interacting system having the same fugacity. The operator 
{}^ then picks up the residues due to interactions. We then may write, 

G (x, x'; r) = {G (x, x'; r)} (/) + G (1) (x, x'; r) , (164) 
where {G}^ can be expanded in terms of {X n }( r K 



B. Trapped virial spectral function up to the second order 



We now calculate the second-order expansion function, which accounts for the leading interaction effect. The next- 
order expansion function, could be treated straightforward using exact three- fermion solutions [73| . The leading term 
of {G-|~|-(x, x'; r)}^ takes the form, 



ze^ {Tr! [ e - w / fc « T e T «# t (x) e - T «*+ (x')] } 



(/) 



(165) 



The trace has to be taken over all the single-particle states (i.e., ijj p with energy e p ) for a spin-down fermion. We 
insert in the bracket an identity J2q \Q) (Q\ = 1; where Q refers to the "paired" state (i.e., $q with energy Eq) for 
two fermions with unlike spins. It is straightforward to show that, at the leading order, 



{G n } 



* T J2{ e ~ ep/kBT+T{£p ~ EQ)F pQ ( x . x ')} (/) 



(166) 



where F p q = J dxidx2'0p(xi) ( I > Q (x, xi) <I>q (x', X2) ip p {x.2)- Accordingly, the leading interaction correction to the 
spectral function, {A (k, ui)}^, is given by, 



z 1 + e 



-hui/k B T 



„ n v 



e p ~ Eq + n) e 



-t p /k B T 



F, 



(I) 



(167) 



where F pQ (k) = f dxdxie 4k r -0p(xi)$ Q (x, x x ). 

In an isotropic harmonic trap with frequency ujt, we can solve exactly the two- fermion problem for relative wave- 
functions and obtain {A(k, w)}^ 7 ' using the above procedure. In the end, we calculate 



2tt 2 



(168) 
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as measured experimentally [4^,|43|- Here, fp (u) = 1/ (e hui '/ fesT + 1) is the Fermi distribution function and 

A^ = 4%/27r/(m 3/2 ^) (w + M - e k ) 1/2 , (169) 

is the spectral function of an ideal, non-interacting Fermi gas. To account for the experimental resolution, we may 
further convolute I(k, u) with a gaussian broadening curve. 
In the BEC limit, we may show analytically that, 



{A} {1 > fp ex exp 



-P (V £ k - w - M + E B - \/ek 



(170) 



where Eb — —h 2 /(ma 2 ) is the binding energy. Thus, at large momentum k the intensity due to interactions peaks at 
to + [i = — £k + Eb, with a width ~ y/ksTe^. At low temperatures, the width should be replaced by y/Epe^, where 
the Fermi energy i?p provides a cut-off to the thermal energy ksT. 

We may also calculate the momentum distribution /Oa-(k) = f_°° dojA(k, uj)fp(uj). At large momentum, we confirm 
the Tan relation [f37|], Per(k) ~ X/fc 4 , where the contact I is given by, 



4irz d 



k B T 
far 



(0) 



(171) 



Here, <f) re i is the relative radial wavefunction of the paired state with energy e re i [73| . At low temperatures, a finite 
contact therefore implies a finite spectral weight below the chemical potential. 

In the calculation, consistent with the leading order expansion in A (Is., w), we determine the fugacity from the 
number equation N = — d[Ail + fl^]/d/j,, by expanding the interacting part of thermodynamic potential, Af2, up to 
the second-order virial coefficient. 
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Figure 26: (color online) Contour plots of the occupied spectral intensity at crossover. The intensity = 
A(k, ui) fF(ijj)k 2 /(27r 2 ) increases from blue (10~ 3 I ma x) to red (I ma x) in a logarithmic scale. The calculations were performed 
with harmonic traps at T — 0.7Tf and l/(kFa s ) = +1, 0, —1, with a resulting fugacity at the trap center of z ~ 0.14, 0.42, 
and 0.48, respectively. From ref. [z3|; copyright (2010) by APS. 



Fig. l26l shows contour plots of the occupied spectral intensity of a trapped Fermi gas in the crossover at T = 0.7Tp. 
At this temperature, our results are quantitatively reliable. We observe that, in addition to the response from coherent 
Landau quasiparticles (black lines), there is a broad incoherent spectral weight centered about = — e^+Es (white 
dashed lines), where et = h 2 k 2 /(2m) and Eb = —h 2 /(ma 2 ) is the binding energy. Thus, the spectra clearly exhibit a 
gap-like double peak structure in the normal state. This is a remarkable feature: the dispersion at negative energies 
seems to follow the BCS-like dispersion curve, u> = — sj (e^ — /i) 2 + A 2 , and behaves as if the gas was superconducting, 
even though we are above the critical temperature T c . Therefore, the incoherent spectral weight indicates the tendency 
of pseudogap: the precursor of fermionic pairing due to strong attractions, i.e., it arises from the atoms in the paired 
state or "molecules". The pairing response is very broad in energy and bends down towards lower energy for increasing 
k. At large l/(kpa s ), the width is of order ^/maxl^T, Ep}ei i . The incoherent spectral weight found by our leading 
cluster expansion is a universal feature of interac ting Fermi gases. At large momentum k kp, it is related to the 
universal 1/fc 4 tail of momentum distribution (67L Il44j |. 
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C. Comparison of theory with the JILA experiment 
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Figure 27: (color online) Single-particle excitation spectra on the BEC side of crossover, a-c, Cluster expansion predictions 
(z ~ 0.1 and fi ~ — 1.08-Ef). d-e, Corresponding experimental data [4tj |. a, The linear-scale intensity map. Our results were 
convoluted with a gaussian broadening curve of width a = 0.22Ef, to account for the measurement resolution [46| . The black 
line shows upper free-atom dispersion. The red dashed line is the lower dispersion curve of molecules, obtained via fitting 
each fixed-fc energy distribution curve (in b) with a two gaussian distribution. It agree fairly well with the experimental result 
(white symbols), b, Energy distribution curves for selected values of k. c, The occupied density of state (DOS). Blue dashed 
lines show the experimental peak positions. Adapted from ref. [76| ; copyright (2010) by APS. 

For a close comparison with experiment j46| . we perform calculations using realistic experimental parameters, 
including the measurement resolution. Fig. [2T]presents the results on the BEC side of crossover with l/(kpa s ) = 1.1. 
The temperature T = 0A5Tp is estimated from an initial temperature Tj = 0.17Tp obtained before the field sweep to 
the BEC side [46[. The experimentally observed upper and lower features, caused respectively by unpaired atoms and 
molecules, are faithfully reproduced. In particular, the experimental data for the quasiparticle dispersion of molecules, 
marked by white symbols, agrees with our theory (lower red dashed line). There is also a qualitative agreement for 
the energy distribution curves (Figs. l2~Tb and \Tfb) and the occupied density of states (Figs. l2~Tb and |2~TD . A narrow 
peak due to free atoms and a broader feature due to molecules are reproduced theoretically with very similar width 
at nearly the same position. It is impressive that the simple quantum cluster expansion is able to capture the main 
feature of the experimental spectra. 
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Figure 28: (color online) Single-particle excitation spectra of a strongly interacting Fermi gas. a-c, Cluster expansion predictions 
(z ~ 6 and fi ~ 0.37 Ef)- d-e, Corresponding experimental data [4q |. In e, for the experimental energy distribution curves, we 
use a larger value of k (i.e., enlarged by a factor of 5/3) to account for a scaling discrepancy due to many-body correlations. 
Adapted from ref. 0; copyright (2010) by APS. 



Fig. [25]reports the spectra in the unitarity limit at the critical temperature T c ~ 0.22V 
the use of a cluster expansion becomes highly questionable as the fugacity at the center z 



At such low temperatures, 
^ 6 >> 1. Nevertheless, we 
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find that the dispersion curve is lowered by the attractions by an amount comparable to the Fermi energy ep , as shown 
clearly by the red dashed line in Fig. l2"5k . The calculated energy distribution curves bifurcates from a single peak 
with increasing k and becomes dominated by the lower molecular branch (Fig. l28b). which eventually leads to the 
bending back of the dispersion curve to negative energy. This picture may be view as an indication of the existence of 
a pseudogap, which is consistent with the experimental findings (Fig. 128b ). This surprisingly good agreement merits 
further investigation. We conjecture that even at these relatively low temperatures the virial expansion captures the 
dominant two-body correlations measured in these experiments, apart from a possible overall scaling factor due to 
the missing higher-order terms. 



VI. VIRIAL EXPANSION FUNCTION AND WILSON COEFFICIENT 



In this section, we discuss briefly the relation between virial expansion and Tan relations, both of which provides 
useful insights to the challenging many-body problem. The virial expansion is a natural tool to bridge few-body 
and many-body physics, while the exact Tan relations give perspective from the point of view of short-distance 
and/or short-time scale. It has been shown by Braaten and Platter that Tan's rela tions can be understood using the 
short-distance and/or short-time operator product expansion (OPE) method (69l Il46l |. in which the few-body and 
many-body scales are separated. At this point, there should be a close relation between virial expansion and Tan 
relations. Here, we show that the Wilson coefficient appearing in the OPE equations is given by the virial expansion 
function (l2| . 



A. Operator product expansion method 



The OPE gives a powerful tool to understand the strongly correlated many-body system in the short -dist ance/short 



time limit. It is a hypothesis independently conjectured by Wilson, Kadanoff, and Polyakov in 1969 [145| . The OPE 
expands the product of local operators at different space-time points in local operators with coefficients that are 
functions of the separation in space and time. For density correlation, it takes the form, 

p a (x, r) p a , (x', rO = W ™> ( x - x '> * ~ t') °c, (172) 

c 

where the sum is over infinitely many local operators Oc [(x + x Q/2 , (t + rO/2] and W^ a , (r — r', r) are called Wilson 
coefficients. The original hypothesis concerns the real time t [145] . Here, we generalize it to an imaginary time r 
via the analytical continuation, t = —it. As a result, the Wilson coefficients defined in this way are amenable for 
calculations at both zero and finite temperatures. The Wilson coefficients rely only on few-body physics. Hence, in 
order to determine W^ a , of a local operator Oc at zero temperature, one may choose a simple few-body state for 
which (Oc) 7^ and match the expectation values on both sides of Eq. (j!72[) . At finite temperatures, however, this 
matching procedure may be considerably complicated. 

In the short-distance/short-time limit, only a few t erms in t he sum of Eq. (|172|) contribute. By neglecting the 
un- important single-particle contribution, it was shown |l46l . ll47l | that after a Fourier transform (q — > oo and uj — > oo), 



S aa . (q, w, T) - (q, u>, T) ~ W aa - (q, u, T) X, (173) 

where X is the Tan' s contact. At zero temperature, the Wilson coefficient of the DSF has been determined by Son 
and Thompson |l47| . 



B. Wilson coefficient from the virial expansion function 



The relation between the virial expansion function and the Wilson coefficient becomes evident, if we expand both 
sides of Eq. (|173[) in fugacity. As W aa i involves only the few-body physics and hence does not contain the fugacity z, 
a count of the term z n on both sides of Eq. (|173p leads to 

W aa , (q,w,T) = ^AS ro ,, 2 (q, W ,T) (174) 
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and 



AS m ,. n (q, u, T) = —AS aa ,. 2 (q, u, T) , 

C2 



(175) 



where I2 = z 2 lQn 2 V C2 / is the contact up to the second order expansion. Therefore, the Wilson coefficient is given 
by the second expansion function, in the case of two-body contact interactions. This result is obtained by applying the 
OPE and virial expansion method. As a result, in principle it should be valid at temperatures above the superfluid 
transition. However, we may expect that it holds at all temperatures, as both the Wilson coefficient and second 
expansion function are irrelevant to the many-body pairing in the superfluid phase. The many-body effect enters 
through the many-body parameter of contact only. As shown by Eq. (|175|1 . in the limits of q — > 00 and 10 — > 00, the 
virial expansion functions becomes proportional to the contact coefficients, as a direct result of the OPE hypothesis. 
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Figure 29: (color online) Wilson coefficients f aa i = V mhw z ^ 2 (z 2 /T2) AS atT ' t 2 at g = 3, 5 , and 10. With increasin g m omentum 
and/or frequency, / CT(T / approaches smoothly to the T — result by Son and Thompson |l47l |. Adapted from ref. [13 |. 



At zero temperature, the Wilson coefficient of the DSF of a unitary Fermi 
using the matching procedure using diagrammatic theory. It is given by [1 

f n /(Vmf^ 3/2 ), where, 



as ca n be analytically calculated, by 
147j WTf° = f n /(V^huj^ 2 ) and 



/ft 



1 y/l-x/2 1 1 

4tt 2 (l _ x f 4tt 2 2x1/! - x/2 



m 2 l 



l-x 



ir 2 <d{x-l) 



(176) 



and 



■hi - 



1 1 

4vr 2 J2x~ 



In 



\/2x — x 2 



4tt 2 2x^/1 - x/2 



In' 



1 + V2x 



1 



-TT 2 0(X-1) 



(177) 



x = h i 2 q 2 1 (2mtkjj) , and O is the step function. On the other hand, the second virial expansion function of the DSF 
at zero temperature can be calculated from the trapped results in the limit of large q = [h 2 c[ 2 / (2mfcsT)] 1 / 2 . 

In Fig. [2Uwe check the validity of Eq. (|174|) at zero temperature, by calculating (mfi,) 1 / 2 w 3 / 2 (z 2 /Z2)A5 crcr ' ! 2 at 
different momenta. With decreasing temperature T or increasing q, it approaches gradually to (mh^^uj^^Wj^J when 
Co > 5 /2. This confirms numerically that Eq. (|174p holds at zero temperature at large momentum and frequency. 
For small frequency (i.e., Co — > 0), the Wilson coefficient becomes divergent. The confirmation of equivalence in this 
limit is stringet and requires a large value of q. Our virial expansion function at q up to 10 is unable to approach 

and W^ =0 have an 



the Wilson coefficient at Co < q /2. We also note that, in the limit of large frequency, W, 
interesting high-frequency power-law tail io~ b / 2 |l47l Il48l |. 
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Ab n (theory) 


Ab n (experiment) 


c n (theory) 


2 


1/^2 [97] 




1/tt [77>i29] 


3 


+1.05 ±0.01° 
-0.35510298 6 
-0.3551030264897° 
-0.3573 ± 0.0005 d 
-0.3551 ±0.0001 e 


-0.35 ±0.02 [30J 


-0.1408 ±0.0010 7 
-0.1399 ± 0.0001 9 


4 


-0.016 ±0.004 [78] 


0.096 ± 0.015 h 
0.096 ±0.010* 




.1 


0.0017 < Ab 5 < 0.101 [78] 







Table I: List of the theoretical predictions and experimental measurements for the virial coefficients of a homogeneous Fermi 
gas i n the unitary limit. T he la st co lum n shows the c onta ct coefficients. For the superscripts (a)-(i), the references are: (a) 

ma, (b) q, (c) a, (d) ma, (e) ma, & mi, ( g ) ma, w a ^ $ m. 

This is fairly evident in the second order virial expansion functions. 

The identification of the Wilson coefficient as the virial expansion function is very useful. For example, in the system 
where the three-body interactions dominate, we anticipate that the third virial expansion function would give the 
Wilson coefficient. At this point, we note that, for identical bosons with a large scattering length in which three-bod y 
Efimov physics occurs, the Wilson coefficient and new universal relation have been derived very recently |l49l . Il50l |. 

VII. OUTLOOK 

In this review, we have demonstrated that virial expansion provides a powerful tool to understand a normal, strongly 
correlated atomic Fermi gas at temperature down to a half of the Fermi degenerate temperature. The virial predictions 
generally agree well with the experimental measurements. 

A. Successes 

These remarkable results cover both static and dynamic properties. 

• For thermodynamics, in the calculation of virial coefficients in the strongly-interacting regime, a convenient 
way is proposed, based on the few-particle solutions in harmonic traps. The exact three-fermion solu- 
tion leads to a very accurate determination of the long-sought third virial coefficient in the unitary limit: 
A63 = —0.3551030264897. The resulting virial equation of state serves as an important benchmark for accurate 
experimental measurements (see Sec. II). It also provides a possible thermometry for strongly interacting Fermi 
gases [32]. The calculation of the fourth virial coefficient in the unitary limit has been attempted, by solving 
numerically the four-fermion problem. In Table 1, we summarize the past theoretical and experimental efforts 
in determining the virial coefficients of a strongly interacting homogeneous Fermi gas. 

• For the universal Tan's contact I that governs the short-range/short-time physics, thanks to the adiabatic 
relation, we can virial expand it in terms of contact coefficients. In the unitary limit, the second and third 
contact coefficients, Ac 2 = 1/tt and AC3 = — 0.1399 ± 0.0001 (see Table 1), provide a good explanation for the 
recent measurement in harmonic traps (Sec. III). 

• For dynamic properties, the dynamic structure factor and single-particle spectral function can be virial expanded 
as well, in terms of virial expansion functions. The second order expansion (in the leading order of interactions) 
gives a good qualitative understanding of recent experimental measurements on two-photon Bragg spectroscopy 
(Sec. IV) and momentum-resolved rf-spectroscopy (Sec. V), at temperature down to the onset of superfiuid 
transition. 
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B. Future developments 

Encouraged by these remarkable achievements, we may foresee a number of potential developments and applications 
of virial expansion in the near future. 

1. Higher-order expansions and new applications 

It is technically straightforward to calculate higher-order virial coefficients and expansion functions. However, 
much heavier numerical efforts would be involved. Owing to the ever-growing power in computation, we anticipate 
optimistically that the fourth and fifth virial coefficients could be calculated accurately. Accordingly, the third to fifth 
virial expansion functions for the single-particle spectral function may be determined. These results will clearly bring 
in-depth understanding of the existing measurements on thermodynamics and spectral function. 

In the novel atomic systems such as multi-component Fermi gases or strongly interacting Bose gases, where the 
three-body or four-body physics becomes important, virial expansion would be particularly useful. Using the third or 
fourth expansion functions, we anticipate to address the many-body consequence of the multi-component (i.e., triplet) 
pairing and Efimov physics. New universal relations may be predicted. 

With these in-mind, we note that the methodology of virial expansion is very general. It can be used as well to study 
many other interesting properties of strongly-correlated atomic Fermi gases, which now become available with current 
experimental techniques. Important examples includes the universal transport coefficient (i.e. the shear viscosity) 
of a unitary Fermi gas, which has been investigated already by the damping rate in collective excitations and by 
the hydrodyn amic expansion (33| . fermionic pair ing i n low - dimensions, which can be probed by the rf-spectroscopy 
|38l . 1 15 lL Il52| , and non-s-wave fermionic pairing 

EE! Em. 

2. Insights for reliable low-temperature strong -coupling theories 

An important motivation of the virial expansion study is to gain insights for developing reliable low-temperature 
strong-coupling theories. Ideally, we wish to apply in a quantitative manner the virial expansion down to the superfluid 
transition temperature. However, in the deep quantum degenerate regime, where the fugacity is larger than unity, we 
may not anticipate convergence of the virial series, evaluated up to certain order. To extract the infinite-order result, 
it is necessary to apply some resummation techniques [64j . 




z = exp[u/(£ B 7)] 

Figure 30: (color online) Universal ft-function as a function of fugacity. The Pade [2/2] approximant is compared with the two 
experimental data sets, from Salomon's group at ENS (empty squares) (3p| and from Zwierlein group at MIT (solid circles) 
[32i |. The vertical shaded line shows the critical fugacity for the onset of superfluid transition. 

As an interesting example, here we discuss briefly the Pade resummation method, in which a virial series, i.e., the 
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Figure 31: (color online) Universal equation of state obtained by using the third-order virial expansion within the [1/1] Pade 
approximant. The result is contrasted with the experimental data from MIT [3^ |. Here Eq = (3/5)NEf is the ground state 
energy of an ideal, non-interacting Fermi gas. The vertical shaded line indicates the critical temperature for the onset of 
superfiuid transition. 



universal ft,- function Eq. (|97|) . is written into the form, 



,(*) = 



Po + piz + p 2 z z H h p m z"' 

1 + q\z + q-2z 2 H 1- q n z n 



(179) 



This is the so-called Pade approximant of order [m/n]. To the order [1/1], the three Pade coefficients p , pi, and q\ 
can be uniquely determined using the three virial coefficients, i.e., 



hW(z) = 

Pade \ ' 



1 + 


[# 


] + Ab 2 - A63/A62 


z 


1 + 


b ( 2 1] - Ab 3 /Ab 2 


z 



(180) 



In Fig. 1301 we compare the universal /i-function in the Pade [1/1] form with the experimental data. It agrees very 
well with the latest measurement (the MIT data set) reported by Ziwerlein's group at all temperatures. The relative 
discrepancy is about 10% in maximum, comparable with the discrepancy of the two experimental data sets. In Fig. 
1311 we compare the virial equation of state, calculated using h£/}\ with the MIT data set. The third-order virial 
expansion within Pade approximant works extremely well, for temperatures down to the onset of superfiuid phase 
transition, T c ~ 0.16Tp. This remarkable agreement, over a wide parameter window in fugacity, is entirely unexpected, 
since the Pade approximant is not controllable and therefore its application can not be justified a prior. We anticipate 
that the accuracy of the virial equation of state could be improved by the inclusion of more Pade terms in h r ^ de (z), 
such as the z 2 term. This is straightforward once the fourth and fifth virial coefficients are accurately calculated. 

It is reasonable to anticipate that the similar Pade approximant may work for the single-particle spectral function. 
In this respect, the virial spectral function within the Pade [1/1] or [2/2] approximant could be useful to clarify the 
delicate pseudogap puzzle in a unitary Fermi gas. 



Acknowledgments 

We have benefited from discussions and collaborations with many physicists: here we would like to especially thank 
Hui Hu, Peter D. Drummond, Peter Hannaford, Chris J. Vale, and Eva D. Kuhnle for valuable interactions in recent 
years, and Tin-L un Ho for his continuous encouragement. We also thank Xavier Leyronas for sending his data file of 
A63 in ref. (lip} , Sylyain Nascimbene and Christophe Salomon for providing us the experimental data of the universal 
/i-function in ref. [3fJ, and Martin W. Zwierlein and Mark J.-H. Ku for providing us the experimental data in ref. 
(32l | . This research was supported by the Australian Research Council Discovery Project (Grant No. DP0984637) 
and NFRP-China (Grant No. 2011CB921502). 



50 



Appendix A: Calculation of C nn i 

In this appendix, we outline the details of how to construct the matrix element C nn ' in Eq. (|T2|) , which is given by, 



Cnn> = J p 2 dpR„l (ft) Rn'l (~) ^2b PI Vl,n>), 



(Al) 



where 



is the radial wave function of an isotropic 3D harmonic oscillator and the two-body relative wave function is 

V 2 f = T{-v hn ,)U{-v l>n ,, |, ^ 2 )cxp(-^ 2 ). (A3) 

Here, for convenience we have set d = 1 as the unit of length. Ln +1 ^ 2 ^ is the generalized Laguerre polynomial 
and U is the second Kummer confluent hypergeometric function. A direct integration for C nn i is difficult, since the 
second Kummer function has a singularity at the origin. The need to integrate for different values of v\ <n i also causes 
additional complications. 

It turns out that a better strategy for the numerical calculations is to write, 



by using the exact identity, 



Therefore, we find that 



where 



k — vi n ' V 2fc! 

k=0 



r(-M-4* 2 ) = £^^. (A5) 

k=0 



k — vi n ' V 2fc! 



r V 1 , r(fc + 3/2) . 



<w* = J ' P 2 d P Rni ( P ) Run (0 R k0 (?yp\ (A7) 

^ ' 

can be calculated to high accuracy with an appropriate integration algorithm. In checking convergence of the sum- 
mation over k, we find numerically that for a cut-off n max (i.e., n,n' < n max ), C l nn , k vanishes for a sufficient large 
k ^> fcmax ^° 4yz max . 

In practical calculations, we tabulate C l nn , k for a given total relative angular momentum. The calculation of C nn > 
for different values of i>i >n i then reduces to a simple summation over k, which is very efficient. Numerically, we have 
confirmed that the matrix C nn i is symmetric, i.e., C nn i = C n / n . 



Appendix B: Calculation of Sj >n 

The calculation of s; „ seems straightforward by using the Bethe-Peierls boundary condition in hyperspherical 
coordinates (|77|) . However, we find that numerical accuracy is low for large n and I due to the difficulty of calculating 
the hypergeometric function 2 -F\ accurately using IEEE standard precision arithmetic. We have therefore utilized 
MATHEMATICA software that can perform analytical calculations with unlimited accuracy. For this purpose, we 
introduce As;.„ = s;. n — si tn . After some algebra, we find the following boundary condition for t = As; jn /2, 

T{-l) n+l T{n + l + l + t) 



sin(Trf) = W-V^? — ' —rf(t), (Bl) 

l/ 3 2T(l + f)r(n+l + t) 
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where we have defined a function 

/ (t) = 2F1 (-n - t,n + 1 + 1 + 1, 1 + |; ^ . (B2) 

The above equation can be solved using the MATHEMATICA routine "FindRoot", by seeking a solution around t = 0. 
It is also easy to write a short program to solve Eq. (|B1[) continuously for 71 < n max = 512 and I < l max = 512. In a 
typical current PC, this takes several days. The results can be tabulated and stored in a file for further use. 
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